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SUMMARY 


The effect of steady state circumferential total pressure distortion on 
the loss in compressor stall pressure ratio has been established by analytical 
techniques. Full scale engine and compressor/fan component test data were 
used to provide direct evaluation of the analysis. Favorable results of the 
comparison are considered verification of the fundamental hypothesis of this 
study. Specifically, since a circumferential total pressure distortion in 
an inlet system will result in unsteady flow in the coordinate system of the 
rotor blades, an analysis of this type distortion must be performed from an 
unsteady aerodynamic point of view. By application of the fundamental 
aerothermodynamic laws to the inlet/compressor system, parameters important 
in the design of such a system for compatible operation have been identified. 

A time constant, directly related to the compressor rotor chord, was found to 
be significant, indicating compressor sensitivity to circumferential dis- 
tortion is directly dependent on the rotor chord. 

As an initial step in the investigation of the effects of time dependent 
total pressure distortion on the compressor stability characteristics, an 
analytical model of turbulent flow typical of that found in aircraft inlets 
has also been developed. Due to the non-deterministic (random) nature of 
this type of flow distortion, the flow analysis requires use of statistical 
methods. These methods were combined with basic fluid dynamic concepts 
to provide a usable analysis technique. With this model, the power spectral 
density function and root mean square level of the time dependent total 
pressure take on considerable significance as indicators of the strength and 
extent of low pressure regions that are important in the compressor reaction 
to inlet flow disturbances. Spectra obtained from the model were compared 
with those obtained in tests of a Mach 3 mixed compression inlet to illustrate 
the technique of determining the mean size and strength of instantaneous low 
pressure regions by statistical techniques and to verify the turbulent flow 
model. Excellent agreement was obtained in the comparison verifying this 
fundamental approach. 

Both the steady state distortion/compressor analysis and the turbulent 
flow model are considered developed to the point necessary to initiate the 
development program to achieve the long term program objective of combining 
these results to establish a fundamental relationship between both inlet 
steady state circumferential distortion and turbulence and loss in compressor 
stall margin. 
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INTRODUCTION 


i Preceding page blank 


Inlet /engine system stability problems have grown to major proportions 
with the continuing press to improve performance and reduce system weight and 
volume. The need to solve such problems and to understand the effect of 
inlet total pressure distortion on engine compressor stability has become 
critical. To date, solutions to the problem of inlet/engine compatibility, 
have had to come from experimental results since adequate stability analysis 
methods were not available. This has resulted in extensive inlet and engine 
test requirements. Notwithstanding, the important design variables for 
inlet/engine stability remained obscure. 

An analytical approach that considers the fundamentals of the dynamic 
interaction between inlet flow and engine compressor is needed to augment the 
use of the traditional empirical distortion factors. The method needs to 
be sufficiently detailed to provide insight into the basic interaction and 
yield workable accuracy, yet not detailed to the point of being expensive 
and cumbersome to apply. 

This program, initiated in April 197? , has been oriented toward 
developing basic relationships between inlet flow distortion and turbulence 
and the loss in compressor stall margin. A five task approach has been 
established. The initial two phases, which comprise the subject matter of 
this report, were designed to develop the fundamental techniques required 
for successful completion of the program. Future studies combine these 
fundamental analyses to relate inlet flow distortion and turbulence to the 
loss in compressor stall margin. These analyses can then be used with 
data from existing inlet/engine tests to establish procedures capable of 
predicting compressor stability margin during the design phase of a 
propulsion system. 

The objective of Task I is to develop an analytical technique to relate 
inlet circumferential total pressure distortion to the loss in compressor 
stall margin. A steady state circumferential distortion appears as time 
variant in the rotor coordinate system. The developed analysis is unique 
since it considers the effects of this unsteady flow on the compressor stage 
characteristics. Secondly, the effects of flow distortion are established by 
consideration of only the stall margin changes caused by distortion, 
eliminating need for detailed construction of individual stage and compressor 
performance maps. Favorable comparison between results of the analysis and 
experimental data are considered to have verified this approach. 

The objective of Task II is to develop a statistical model of inlet 
turbulent flow. This was accomplished by the combination of two engineer- 
ing disciplines: fluid mechanics and statistical mathematics. Based on 

the fundamental hypothesis that the time dependent total pressure fluctua- 
tions are a direct result of streamline curvature rather than acoustic waves, 
it was assumed that these pressure fluctuations could be described by a 
random distribution of descrete vortices transported by the mean flow. The 
laws of fluid mechanics were used to describe the fluid dynamic character- 
istics of the vortices, while the statistical methods were used to handle 
the random properties of the flow. Results of the analysis were verified by 
test data. Through this model easily measured inlet flow properties such as 
total pressure RMS level and power spectral density function can be inter- 
preted in a context meaningful to engine stability. 


PRECEDING PAGE BLANK NOT FILMED 



SYMBOLS 
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English 

A 

a 

B 

C 

C 

c 

C 1 

c 2 

D 


d 


d( ) 

E 

ERF 

ERFC 

e 

F 

f 

f( ) 
Qx(f ) 


H 

h 

h 

i 

k 

k 

L 

M 

m 

m 

N 

N 

n 


n 


n 

P 

q D 

R 

RMS 

r 

S 


area 

vortex core radius 

vortex strength coefficient 

coefficient 

constant 

chord 

spring constant 
damping factor 
coefficient in 
transformed hyper- 
beta function 
coefficient in 
transformed 
beta function 
differential of ( ) 
energy 

error function 

complementary error function 

base of natural logarithm 

force 

frequency 

function of ( ) 

real one-sided power 

spectral density function 

inlet duct height 

enthalpy 

probe location 

/-I 

beta density coefficient 
reduced frequency = wc /2 u 
lift 


ASM 

Sx(f) 

T 

t 

U 

U 


u 

V 




( )2 


Greek 

a. 

a. 


Mach number r 

exponent in beta function y 

mass A 

frequency of occurence e 

rotor RPM tj 

direction of vortex 17 

rotation ( + ,-) © 

exponent in beta and hyper- v 
beta density function 
harmonic number 7 r 

pressure 1 2 ^ 

dynamic pressure = -*-p U <7 

ratio * ° r 

root mean square T 


LaPlacian variable 


= loss in stall margin 
= complex power spectral 
density function 
= transfer function 
= time 

= axial velocity (vortex 
analysis) 

= relative velocity (compressor 
analysis ) 

- perturbation velocity in 
x direction 

= vertical velocity (vortex 
analysis ) 

= perturbation velocity in 
y direction 
= resultant velocity 

= coordinates fixed to inlet 
probe 

= coordinates fixed to vortex 

= mean value of ( ) 

= square of mean value of ( ) 

= mean square value of ( ) 


= angle of attack 
= flow angle 
= circulation 

= ratio of specific heats = 1.4 
= difference 

= small distance from probe 
= efficiency 

= total pressure recovery 
= circumferential angle 
= coefficient of kinematic 
viscosity 

= 3.14159 


= density 
= RMS value 
= delay time 

= nondimens ional time = tU/c 
= nondimens ional time constant = 
time constant ^u/c 
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S2 

co 


SYMBOI5 (Continued) 


= vorticity 

= angular frequency = 2 7r f 

English Script 

L .[ ] = Laplaclan operator 

p( ) = probability density 

function of ( ) 

R x (t) = autocorrelation function 
of x 


Subscript 

a = airflow 

a = core size 

avg = average 

ax = axial 

c = circulatory 

c = compressor 

eff = effective 

g = general 

inst » instantaneous 

L = low pressure region 

max = maximum 

N.C. = non-circulatory 

min a minimum 

a a freestream (uniform conditions) 

p = pressure 

RTR = rotor 

r a radial 

T a total pressure 

v = vortex strength 

9 = tangential (circumferential) 

direction 

o a steady-state 
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TASK I 

EFFECT OF STEADY STATE TOTAL PRESSURE DISTORTION ON COMPRESSOR STALL MARGIN 


The objective of Task I is to relate inlet circumferential steady state 
total pressure distortion to loss in engine compressor stall margin. An ana- 
lytical technique based on the fundamental aero-thermodynamic laws governing 
fluid flow and engine compressor operation has been developed. The general 
approach is outlined below and the details presented in subsequent sections. 

Distorted inlet flow is composed of total pressure levels both above and 
below the average. These regions correspond to deviations in axial flow 
velocity from the mean. In the rotating coordinate system of the rotor, these 
deviations appear as fluctuations in the stream angle or angle of attack rela- 
tive to the rotor blades. Therefore, the flow over the rotor blades is 
basically unsteady and hence steady state distortion, as well as -unsteady, must 
be analyzed by -unsteady aerodynamic techniques. Accordingly, as a basis for 
the study, the effects of a time varying angle of attack on the lifting 
characteristics of an isolated airfoil are established. The results are then 
applied to a compressor rotor blade and by relating the work done by the 
rotor to the lifting characteristics of the blades, the loss in compressor 
stall margin due to an arbitrary circumferential distortion pattern is 
established. 


Isolated Airfoil Analysis 

The primary objective of this specific item is to establish the effect 
of unsteady airflow on the lifting characteristics, and in particular on the 
maximum lift coefficient, of an isolated airfoil. This will include resolu- 
tion of the effects for arbitrary transients in angle of attack. To accomplish 
this objective, it is first necessary to understand the flow phenomena 
involved in delaying the stall of an airfoil beyond its steady state charac- 
teristics when the airfoil is subjected to an unsteady angle of attack 
and then develop a mathematical representation of the process which can be 
solved for arbitrary, time dependent, angles of attack. 

Effect of Unsteady Flow on Lift . - Lift on an airfoil is a consequence 
of unequal pressures acting on the upper and lower surfaces. In potential 
flow these pressures can be computed from the velocity field by use of the 
equations of motion. In the case of unsteady flow, the lift is dependent 
not only on the instantaneous angle of attack but also on the following two 
factors: (l) the inertia or acceleration of the mass of air in prox imi ty 

of the airfoil and, (2) the shedding of the trailing edge vortex which acts 
as a dissipative force. The phenomena are analogous to the forces and 
acceleration of a damped mass/spring system which can be described by a 
linear second order differential equation. Similarly, the lift of an air- 
foil subjected to an unsteady flow can be described in the same manner. 

As an example, the lift per unit span due to an airfoil undergoing vertical 
oscillations at an angular frequency of is : 

L(t) - £irp^ J + |- UirpcC(k) + |pU 2 cir j a Q 


( 1 ) 
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•where: = virtual mass 

7rpcC(k) = "dissipation constant" 

C(k) - function of reduced frequency, k 
k = 6>c/2U 

Similar expressions govern the response of airfoil lift to a wide variety 
of motions. The unsteady lift equations for the various classes of motions are 
summarized in Appendix A. Airfoil lift characteristics of an oscillating air- 
foil are shown in Figure 1 to illustrate the effects caused by the unsteady 
motion. Analytical results are shown compared with test data from Reference 3. 
The qualitative agreement verifies the classical potential flow analysis. 

The effect of the unsteady motion, illustrated in Figure 1 are directly 
related to the reduced frequency, k, which is an extremely important parameter 
in the analysis of unsteady flow over airfoils. In this parameter the ratio 
of chord to airfoil velocity, c/U, is proportional to the time required for a 
disturbance to pass from the leading edge to the trailing edge of the airfoil. 
The time associated with the disturbance (in this case the oscillations) is 
proportional to l/u> . The reduced frequency, k, can therefore be described as 
the ratio of the time associated with a disturbance (l/w) to the time for the 
airfoil to react to the disturbance. 



Figure 1* Unsteady Lift of Oscillating Airfoil 


Unsteady Flow Model - The response of the airfoil to unsteady motions in 
unstalled flow forms the basis on which to develop the phenomenological model 
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of an isolated airfoil subjected to angle of attack excursions beyond the steady 
state stall limit. This is achieved by modeling the physical mechanisms 
involved with a stalling airfoil via the concept of an effective angle of attack. 

When a airfoil is subjected to unsteady variations in angle of attack, the 
pressure distribution about the airfoil does not correspond to that associated 
with the steady state condition for the instantaneous value of angle of attack. 
This is due to the finite amount of time required for flow about an airfoil 
to adjust to the variations in angle of attack. Flow phenomena requiring 
adjustment include the external flow, shed vorticity, and the boundary layer. 
Initially the flow at the airfoil leading edge experiences the change in angle 
of incidence. At later times this new flow angle is felt at subsequent stations 
along the chord of the airfoil. Therefore an effective angle of attack, a eff, 
is hypothesized which lags the instantaneous angle. This angle accounts for the 
finite time required for airflow adjustment and boundary layer separation to 
occur and is modeled mathematically below to enable prediction of the stalling 
lift coefficient of an airfoil operating in unsteady flow. 

In keeping with the findings of an unstalled airfoil, it is assumed that 
the physical mechanisms are governed by a linear second order differential 
equation. Thus, the relationship between the effective and instantaneous angle 
of attack can be written as: 


d 2 ( a eff a o) . 

d? 

where: 


<7 ♦>> 

1 2 

T = 


^inst = 


a eff = 

a = 
o 


a( ° e ~ d ~ ' ~‘ 1 ' * rr ( ° eet - ° c) ■ t , ( “ lnst - 

non-dimensional time = t(U/c) 

actual (instantaneous) angle of attack at time, t 
effective angle of attack 

angle of attack about which the perturbations occur. 


( 2 ) 


The time constants, t, and t 2 are associated with the airfoil/airflow system 
and are to be established from test results. The equation can be solved by 
LaPlace transform techniques for instantaneous angles of attack that vary as 
simple functions of time. This method of solution and the solutions for a sine 
and ramp change are illustrated in Appendix B. The delay in the effective 
angle of attack resulting from a step increase in angle of attack is shown in 
Figure 2. The dependence of this delay on the respective time constants is 
evident . 
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Figure 2. Response of Effective Angle of Attack to a Step- Change 


The reduction in amplitude and time lag of the effective angle resulting 
from an instantaneous angle of attack having a periodic sine variation with 
time is illustrated in the Figure 3. This is the type pattern that a com- 
pressor rotor blade might experience behind a 180 degree circumferential dis- 
tortion pattern. The ratio of the maximum amplitude of the effective angle to 
the maximum instantaneous angle (Equation 3) is dependent only on the system 
reduced frequency and the two time constants and is designated, f(k). 


(a 


eff 


inst 


- a ) 
o max 

- o ) 
o max 


f (k) 


Since only the ratio of the 


x; 


4k 2 T X 2 ) 


(3) 


a t 


4k 2 T 2 2 ) 


cimum angle is of interest the subscript 
"max" will be dropped. Henceforth, the function f(k) will be understood as 
representing this ratio. The function f(k), shown plotted in Figure k for 
various values of the respective time constants, is used along with a Fourier 
Series to establish the airfoil response to arbitrary variations in the instan- 
taneous angle of attack. 


Extension to Arbitrary Variations in Angle of Attack - To establish the 
rotor airfoil response characteristics to any type of circumferential dis- 
tortion pattern, it is necessary to solve Equation 2 for the effective angle of 
attack given arbitrary variations in the instantaneous angle of attack. This 
will enable the compressor characteristics to be determined as a function of 
the circumferential distortion and the subsequent loss in stall margin esti- 
mated. 


A periodic transient, a ^ nst , can be represented by a sum of sine and 
cosine waves, i.e. Fourier Series. Since the governing differential equation 
(Equation 2) is linear, solutions can be superimposed. Therefore, by repre- 
senting the input transient as a Fourier Series and by the use of superposition, 
a solution for an arbitrary transient can be obtained. 
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The Fourier Series representation is as follows : 

OO 0O 

°inst (0) = “o + XX COS (n0) + XX Sin (n6) (4) 

n = 1 n = 1 


Where: n = the harmonic number 

Of Q = average angle of attack 
a^, b n = Fourier Coefficients 

In practice, the number of harmonics required (n) is determined by the accuracy 
required in approximating the input signal. As an example the Fourier Series 
fit of one cycle of the periodic rectangular pattern is shown in Figure 5(a) 
for 10, 25 and 50 harmonics. 



ANGULAR POSITION, deg 

(a) Fburier Series Fit of Input Angle of Attack 

Figure 5. Fourier Series . Fit of the Instantaneous and Computed Effective Angles of Atteck. 

The effective angle of attack is related to the instantaneous angle for 
each harmonic through Equation 3* If the variation in the instantaneous angle 
of attack has a frequency of f cycles/sec., the angular frequency of oscilla- 
tion, u , is 2 7Tf and corresponds to the first harmonic in the Fourier Series. 
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The second harmonic will he twice 2 7Tf or 47Tf. In general, the angular 
frequency of the nth harmonic will be n(2 7Tf). Equation 3 can now be applied 
to each harmonic as illustrated in Equation 5. 


In general: 


“eff 

“inst 



= f(nk) 


n 


(5) 


The effective angle of attack of the total input signal is found by adding the 
solutions for the individual harmonic as indicated by Equation 6. 


- “0 


where ¥(nk) 


00 


Z f (nk) 
n=l 


a cos 
n 


(n0 + ’J'(nk) + 2 f (nk) b sin (n0 + Y(nk) 

n-1 n 


= tan 


-1 


(2nkx^) + tan 


-1 


(2nkt 2 ) 


(6) 


The results for the rectangular periodic pattern are shown in Figure 5(b) for 
an increasing number of harmonics. Although an accurate fit of the rectangular 
wave requires a large number of harmonics, the effective angle is relatively 
insensitive to this number. 

Airfoil Dyn^c Stall - Stall of an airfoil in unsteady flow occurs at 
higher instantaneous angles of attack than that obtained under steady state 
flow conditions. This is indicated schematically in Figure 6(a), where the 
point "D" represents the instantaneous stall point and "B" the steady state 
stall point. This concept results in a time lag in the airfoil response to the 
unsteady airflow and a reduction in the maximum effective angle of attack. 

Both of these items are due to the finite time required for the airflow about 
the airfoil to adjust. This lag in response is indicated in Figure 6(b) for 
a sinusoidal variation in angle of attack and superimposed on the airfoil 
characteristic in Figure 6(c). The relationship governing this effective angle 
of attack is given by Equation 6. It is hypothesized that when a e ff e( l' ua - 1 - 
to the steady state stall value, stall during unsteady flow. will occur. Thus, 
in Figure 6(b) when o; reaches the steady state stall line (line B) the 
airfoil will stall. This stall condition, a eff = a sss is represented for a 
sinusoidal oscillation by Equation 7. 


a sss _ ~ JL o s f(k) ( 7 ) 

°inst °o 

Solution of Equation 7 for the instantaneous angle of attack will yield the 
maximum allowable value for the specific, f(k). Thus: 


‘■inst - O ■ (a - « ) / f(k) 

Inst o sss o 

max 


( 8 ) 
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(a) Typical Airfoil Lift Characteristic 




(c) Effect of Sinusoidal Oscillations 
on Airfoil Lift 


O - steady state (or mean operating point) 
B - steady state stall 
C - maximum instantaneous excursion 
D * "instantaneous" v stall point 
E • maximum effective angle 


Figure 6. The Effect of Sinusoidal Oscillation on Airfoil Characteristics 


The increase in maximum (stalling) angle of attack of an airfoil will there- 
fore be: 


Acs = a . 
max inst 


max 


°sss “ ^ a sss " a o )[ f?k7 


(9) 


This will be the increase in the stalling value of a ^ . as indicated by- 
point D in Figure 6(a). The function f(k) is dependent on the respective 
system time constants, and and the reduced frequency, k. 
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To establish an estimate of the time constants a limited literature survey 



of the effect of unsteady flow on the maximum lift coefficient of an airfoil 
was conducted and is presented in Appendix C. Results indicate that the time 
constants are approximately equal and on the order of 3-5c/U. 

In summary, it was found that the response of a lifting airfoil to an 
unsteady change in angle of attack was in general governed by a second order 
linear differential equation. To represent this unsteady process which is a 
function of the time required for airflow accelerations, shedding of necessary 
trailing edge vortices, and the delay of boundary layer separation, an 
effective angle of attack was hypothesized. By means of this effective angle 
of attack a mathematical representation of the increase in stalling lift co- 
efficient is established by solution of the governing differential equation. 
This is considered an important development since it enables the response 
characteristics of a rotor airfoil subjected to unsteady flow conditions to 
be determined. These characteristics can then be incorporated into a com- 
pressor analysis. 


Compressor Analysis 

The response of a compressor rotor to circumferential total pressure 
distortion will be established by first relating the change in rotor airfoil 
angle of attack caused by the distortion to the required change in compressor 
pressure ratio. This result will then be combined with the unsteady flow model 
for an isolated airfoil to relate the inlet pressure distortion to loss in 
compressor stall margin. Fundamental to this analysis is the assumption that 
the stage or stages that first cause breakdown or surge in the compressor 
operating in undistorted flow are the same limiting stages causing the compressor 
to stall ■when subjected to a distorted flow. This assumption enables the 
analysis to predict perturbations of the stall line due to distortion rather 
than an absolute stall margin level, which would require a stage by stage 
analysis. 

Relate Distortion to Blade Lift Coefficient and Compressor Work . - The 
object of the following development is to relate the total pressure distortion 
at the compressor face to the required additional compressor pressure ratio 
and rotor blade lift coefficient. This is accomplished by means of the follow- 
ing approach. 

The overall performance of a compressor is represented by a compressor map 
as shown schematically in Figure 7. To minimize weight, the engine is designed 
to operate at high stage loadings, near the stall line as shown. When the 
compressor is subjected to a distorted flow, the average work done by the 
compressor on the airflow remains constant, and corresponds to point 0 in 
Figure 7. However, that section of the compressor operating in the region of 
low inlet total pressure must operate at a higher pressure ratio (point 1 in 
Figure 7) to pump the flow to the uniform compressor exit pressure. The 
opposite condition holds for the high pressure regions, which correspond to 
point 2 in Figure 7. The low pressure regions are of prime interest since they 
tend to reduce the compressor stall margin. The additional work required in 
the low pressure regions is assumed to be evenly divided among the compressor 
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stages. For each stage, the relationship between the change in rotor work 
due to distortion and the change in rotor blade lift coefficient can be obtained 
by equating the change in work done on the air to the change in the rotor lift 
charact eristics. This is developed in detail in Appendix D with the following 
result . 



Stall Line 
Operating Line 


In essence the fractional change in work done by each rotor on the airflow 
equals the fractional change in blade lift coefficient or angle of attack. 
Furthermore, the required work increase can be related to the required increase 
in compressor pressure ratio as indicated by Equation 11. 


d( Ah ) 
~ Ah 


r . - 1 

I dRp 
±=I R p 
1 - Rp V 


where 7 


ratio of specific heats = 1.4 


R = compressor pressure ratio 
P 

Ah = rotor work increase 


( 11 ) 
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By combining equation (10 ) and (ll) the increased blade lift coefficient 
and/or angle of attack is found to be: 


da 

a 



1 


I 



1=1 R t 


V 


( 12 ) 


The increased pressure ratio required of the compressor is the negative 
of the change in inlet total pressure due to distortion, or 

( 


Combining this with Equation 12 produces the desired relationship between the 
change in rotor lift coefficient (dC L ), and angle of attack (da), and the 
inlet flow distortion (dP T p) • 



da 

a 


1=1 

T (1-y)/y 

1-R 

p 



m 


This result can then be combined with the effects of unsteady flow on the 
stalling lift coefficient to establish the effect of distortion on the loss in 
compressor stall margin. 

Relate Inlet Pressure Distortion to loss in Compressor Stall Margin . - 
The procedure to establish the loss in stall margin is developed with the aid 
of Figure 8. The steady state and dynamic rotor airfoil characteristic axe 
shown in Figure 8(a). The dynamic characteristics are typical of that pro- 
duced by a circumferential distortion. The actual or instantaneous angle of 
attack on the rotor and resultant lift coefficient are shown as the outer 
ellipse. The maximum operating point is designated point C. The effective 
angle of attack as defined in the unsteady analysis lags the instantaneous 
angle and is shown as the inner ellipse with a maximum at point E. 
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dC L (y-D/v fV 

C L i.RpU-wr Rp 



Figure 8. Transformation from Airfoil Characteristic to Compressor Characteristic 


Equation 12 can be used to transform these changes in rotor character- 
istics of a single stage to the map characteristics of a multi-stage com- 
pressor. This relationship between the required change in compressor pres- 
sure ratio due to distortion and the change in rotor airfoil angle of attack 
is given by 


da (r-i)/y 

— " (1-yW 

1-Rt, 


= Wpj ttp 
R-D f( ^ j RP 


(15) 


On a finite basis, assuming small changes, the airfoil characteristics can be 
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converted to the compressor characteristics hy use of Equation 15. The results 
are shown in Figure 8(b). The relationship between the effective and instan- 
taneous pressure ratios is given by 


(%eff ” R P0) _ ( a eff " a 0 } 

r^Pinst - Hpo) f(P p ) (« inst -a 0 ) 

The case of interest is where Rp inst is a maximum. 

Figure 8. Thus, Equation l6 becomes: 

E Peff ~ E P0 = E PE ~ E P0 _ a E ~ a 0 (17) 

R Pinst” *70 E PC “ E P0 a c " a 0 

Defining the compressor inlet and exit stations as station 2 and 3» respectively, 
Bpp can be established from the distortion level as follows: 


( 16 ) 

This occurs at point "C", 


E PC " P T3 //P T2min ( ^ 3/^2 ^ ^ ^ ^in^T ^ 

R - Rr 0 

150 ' P - P . 

1 - T Tmxn 


(l8a) 


Defining the magnitude of distortion, D.ist, as(P T - P T /P T ) , Equation 18(a) 
becomes min 

Rt 


R. 


P0 


PC 


( 18 b) 


1 - Dist 

Referring to Figure 8(b), the stall margin (SM) with distortion will be: 


(sm) dist = !ro_Jpg 


R. 


( 19 a) 


K) 


whereas stall margin with zero distortion is : 


(SM) 


CLEAN 


= E PB " E P0 


( 19 b) 


R 


PO 
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Therefore, the loss in stall margin (4SM) due to distortion will be the 
difference of Equations 19(a) and 19(b) or 


R 


ASM = 


PA 


= R PE " R P0 


P-nr-i - Rr 


PB PO PB PE 


R. 


PA 


PO 


(20a) 


For a sinusoidal variation in angle of attack, the ratio of the maximum effec- 
tive to maximum instantaneous angle of attack is equal to f(k) as defined in 
Equation 3. Therefore under these conditions Equation 17 can be written as 

R PE " R P0 = f ( k )( R pc “ R Po^ 

This can be combined with Equation 20(a), resulting in the following: 


A SM = f (k) 


/ R PC ■ R P0 
{ *P0 


The pressure ratio, Rp^, was related to distortion in Equation 18. Incor- 
porating this expression into the above, the loss in stall margin becomes: 


A SM = f (k) 


Asm = f (k) 


R 


PO 


l-Dist 

R P0 

1 

l-Dist 


- R. 


PO 




By a series expansion of j- _ - p - ist , A SM can be written 


(20b) 


A SM = f(k) (Dist + Dist + — ---) 


The second term is small for reasonable values of distortion and the relation- 
ship between the loss in compressor stall margin and distortion becomes sinply: 


A SM 
Dist 



A SM 

P Ttein ) /P T 


= f(k) 


( 21 ) 
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The loss in compressor stall margin can be established by use of Equation 
21 for a 180 degree sinusoidal circumferential distortion pattern. Results of 
this computation are shown in Figure 9 for T ^ = r 2 and f° r various values of 
the non-dimensional time constant* The abscissa has been modified to include 
the number of lobes, n, in the circumferential distortion pattern. The graph 
is thus generalized to enable the stall margin loss to be found for a com- 
pressor of reduced frequency, = (Jc/2U and for single or multiple lobe 
distortion patterns. This loss for several typical distortion patterns is 
shown in Figure 10 for three assumed compressors having reduced frequencies, 
k c , of 0.05, 0.10 and 0.15, respectively, and time constants r ± = = 3.5. 



MODIFIED REDUCED FREQUENCY, k*nk -nw c/(2U) 

c c 


Figure 9 The Effect of Compressor Reduced Frequency 
and System Time Constants on Loss in Normalized 
Compressor Stall Margin 





Figure 10 Tolerance to Sinusoidal Distortion for 
Different Compressors 



The difference in stall margin loss resulting from the assumption of first, 
second, third or fourth order systems (i.e., solutions to a first, second, 
third, and fourth order differential equations) is demonstrated in Figure 11. 
This knowledge will be used later in comparison of analytical predictions with 
actual compressor test results to obtain the proper time constants. 



Figure 11 A Comparison of a First, Second, Third, and Fourth Order System 
on Loss In Compressor Stall Margin 

i 


The loss in compressor stall pressure ratio resulting from arbitrary 
(non-sinusoidal) circumferential distortion patterns is determined by use of 
the Fourier Series techniques. A computer program has been developed to 
mechanize the calculation and is documented in Appendix E. Square and rectan- 
gular wave patterns can be evaluated by this technique. The loss in stall 
margin for a l 80 ° square, wave distortion pattern as opposed to a l 80 ° sine 
wave pattern is shown in Figure 12. The results are shown in terms of a general 




reduced frequency k which is equal to twice the product of the compressor 
reduced frequency, the number of lobes in the distortion pattern, n, and 
the time constant T . This type of presentation is valid only for the case 
where = r 2 = T . Results indicate that the square wave pattern will cause 
a greater loss in stall margin than the sine pattern. However, a rectangular 
pattern with a sharp edged profile cannot be realized with the mixing asso- 
ciated with non-uniform flow. An estimate of the effect of this mixing on the 
expected stall pressure ratio can be obtained by modifying the square profile 
as indicated in the insert of Figure 12 and computing the loss in stall margin 
caused by such a pattern. The effect of such modifications are also shown in 
Figure 12 . 



Figure 12 The Effect of Distortion Profile on Loss in Compressor Stall Margin 


Although the analysis includes certain simplifying assumptions, it is 
based on fundamentals of fluid mechanics and aerodynamics, such as the compres- 
sor flow/work balance and the approximated lag functions which are known to 
characterize the airfoil and hence stage response to unsteady airflow. The 
analysis has shown that the stall margin loss is directly a function of the 
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distortion level ((P^ - PTmin^/^ 1 )’ the shape of the distortion pattern and of 
the compressor rotor reduced frequency, k . 


Application and Generalized Curves 

The technique relating arbitrary circumferential distortion patterns to 
the loss in compressor stall pressure ratio has been established. This analy- 
sis has been computerized. Documentation and instructions for use of this 
program is given in Appendix E. However, many inlet distortion patterns can 
be approximated by standard patterns of sine or rectangular wave shape. 
Furthermore, the loss in stall margin resulting from such distortion patterns 
is basically dependent upon only the extent of the low pressure region, the 
rotor time constant, T , and the rotor reduced frequency, 1^. As a result, the 
loss in stall margin associated with these patterns can be presented as a 
function of these three variables, ©l, t , and 1^. Therefore a set of general- 
ized curves have been compiled that can be readily used to estimate the loss 
in compressor stall margin for these standard patterns. 

The three basic distortion patterns utilized to compile the generalized 
curves are shown in Table III. Applicable definitions and nomenclature are 
presented in Table I as an aid in estimating the loss in stall pressure ratio 
for those patterns defined in Table III. A means to approximate the compressor 
reduced frequency is shown in Table II. 

Table III is composed of two parts. The curves outlined in the first part 
are completely general and can be used with a non-dimensional time constant 
chosen by the user to establish, for example, the effect of different time 
constants on a comparison between the analysis and a specific set of compressor 
test data. The second part pertains to those generalized curves utilizing a 
fixed non-dimensional time constant with a value of T s 3,5, It will be shown 
in the Data/Analysis Comparison that this constant will produce a good match 
between test results and the analysis. The specific curve applicable to a 
given problem will depend on the distortion pattern and available information 
on the non-dimensional time constant. Table III is intended to give the 
required guidance for use, of the specific curves, contained in Figures 13 
through 16 . 


Comparison of Analysis with Test Data 

A limited comparison between results of the analysis with test data was 
conducted to add credence to the analysis developed herein, which is based 
solely on theoretical grounds. A literature search was conducted and the data 
limited to that readily available in published sources. Much of this data 
lacks specific details and the reduced frequency of the respective turbo- 
machinery has been estimated. A description of each test vehicle along with 
the source of information and the reduced frequency is given in Table IV. 

The data comparison is made assuming a single non-dimensional time con- 
stant to be valid for all turbomachinery. In such a case the compressor 



TABLE I 

DEFINITIONS AND NOMENCLATURE 


Distortion 
SM 
A SM 

A SM /(fp m - P )/? ) 
Al T Tmin '' T' 


co 

N 

c 

U 


Dm > _ 


Lip 


'ax 



- P, 


Tmin 



= Stall Margin 


= Loss in SM along a line of constant 
corrected rotor speed 


= Normalized Loss in Stall Margin 

= number of lobes in distortion pattern 
(Refer to Table III) 


= Compressor reduced frequency =coc/2U 
= Rotor angular frequency =2 tt n/6o 
= Rotor RPM 

= First stage rotor chord at the tip 

= Velocity relative to the first stage rotor 
at the tip 

= Generalized reduced frequency = 2 k c n 
= Non-dimensional time constant* 

= Rotor tip diameter 
= Axial velocity into the rotor 


* It will be shown in the Data/Analysis comparison that a system 
having time constants = t 2 = 3.5 is a good approximation for 
all fan and compressor systems. 
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TABLE II 


ESTIMATING COMPRESSOR REDUCED FREQUENCY 


Compressor reduced frequency is defined as: 
k c = co c/2U =7rNc/60U 

With lack of specific data this reduced frequency can be estimated by 
assuming axial inflow into the rotor. The reduced frequency is established 
by use of the velocity diagram, with the following result: 


k c = dt: 


tip 


, .V ax \2 3 /Vax ,b 

1 ' 1/2 <u — ' + H ( u — } 


RTR 


RTR 


High flow, high tip speed compressors will have a ratio of (Vax/u )‘ 

RTR 

on the order of 0.2. The approximate reduced frequency will then be 
as follows : 

k = .9 c/D, . 

c ' tip 


The range of this parameter normally lies between .05 and .15. 


26 



TABLE III 


USE OF GENERALIZED CURVES TO ESTIMATE LOSS IN STALL PRESSURE RATIO 




NUMBER 

LOW PRESSURE 


JISTORTION PATTERN 

DISTORTION PATTERN 

OF 

REGION WIDTH, 

TIME CONSTANTS FIGURE 

DESCRIPTION 

SHAPE 

LOBES ,n 

®L 

T ^ T 2 NUMBER 


.ulti-lobe Sin*' 
ave 



ulti-lobe Square 
ave 



n 


n 


l8o/n 


180/n 


t t 13 


T T 13 


ingle Lobe Square r q 
ave 1 - ^L 


2 ir 


1 


T T 14 


ulti-lobe Sine 
*ve * 


ulti-lobe Square 
ave * 



ingle Lobe Square 
ave * 



a 

2jr 

n 


e 


l 


9 t 


3.5 3.5 16 


* This data incorporates the preliminary estimate of a universal non-dimensional 
time constant. 
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GENERAL REDUCED FREQUENCY, k g = 2k c nT 

Figure 13. The Effect of the General Reduced Frequency on the Loss in Stall Margin 
for Sine and Square Wave Patterns. 

























NORMALIZED LOSS IN STALL MARGIN 


MODIFIED SEDUCED FREQUENCY, k - nk, - n <J q C/(2U rel ) 


Figure 15. The Effect of the Compressor Reduced Frequency on the Loss in Stall 
Margin for Sine and Square Wave Patterns. 

















TABLE IV 

SOURCE OF TEST INFORMATION 


Description of 
Test Vehicle 

Type of 

Distortion- 

Pattern 

Designation 

Source 

Reference No. 

First Stage Rotor 
Reduced Frequency, 
k°tip 

Figure No 

General Electric 
Compressor 
Component Test 

Sine 

Profile 

GE4 Block I 

4 

— ' . 089 * 

17 

General Electric 
Compressor 
Component Test 

Sine 

Profile 

GEU Block II 

5 

.059* 

17 

General Electric 
Turbojet Engine 

Rectangular 

Profile 

J85-13 

6 

.087* 

18 

General Electric/ 
NASA High Speed 
Single Stage 
Research Com- 
pressor 

Rectangular 

Profile 


7 

OJ 

CO 

O 

* 

18 

Pratt and Whitney 
Aircraft Fan-Low 
Compressor Rig 
Test 

Rectangular 

Profile 


8 

. 10 * 

19 

Pratt and Whitney 
Aircraft 3 Stage 
Low Speed Research 
Compressor Rig 

Sine 

Profile 


1 

.0325 

.065 

.130 

20 

Rolls Royce 6 Stage 
Compressor 
Component 

Rectangular 

Profile 


9 

. 067 * 

21 


* Estimated 



NORMALIZED LOSS IN STALL MARGIN 



Distortion Profile 
(a) GEh - Block I (Reference k) 



Distortion Profile 
(b) GEl* -r Block II (Reference 5) 


Figure 17. Comparison of Analysis With General Electric Full Sc*le 

Compressor Test Data 



Distortion Extent ,0, DEG 
(a) J-85 Engine (Reference 6) 



Distortion Extent ,9* DEG 
(b) NASA Single Stage Rotor (Reference J) 


Figure 18. Comparison of Analysis with Test Data 
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NORMALIZED LOSS IN STALL MARGIN 


©»© > 



DISTORTION EXTENT, 0, DEG 


O '® $ 



DISTORTION PROFILE 


Figure 19. Comparison of Analysis with Pratt and Whitney Fan/low 
Compressor Component (Reference 8) 




Figure 20: 
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Comparison of Analysis with Test Results of the Pratt and Whitney 
Aircraft Low Speed 3 Stage Research Compressor (Reference 1) 





sensitivity to distortion is dependent only on the rotor reduced frequency 
as shown earlier (Effect of Inlet Pressure Distortion on Loss in Stall Pressure 
Patio). The comparison is presented in terms of the normalized loss in stall 
pressure ratio in Figures 17 through 21 for the constant non-dimensional time 
constant of 3.5. Without exception, the trends predicted by use of the analy- 
sis are in very good agreement with results of the tests, and in many cases, 
good quantitative agreement is obtained. (See for example Figure 20) 

The non-dimensional time constant of 3.5 was assumed universally appli- 
cable in the above comparison. As shown it produces reasonable agreement with 
data and as a result is recommended for use during the development phase of an 
inlet/engine program. However, the analysis also provides the mechanism to 
improve data/analysis comparison once test data of production type hardware 
is available. For example, the prediction shown in Figure 21 for T = 3.5 can 
be modified by assuming a time constant, r = 2.5, bringing the analysis in 
exceptionally close agreement with the test results as shown in Figure 22. 

This further demonstrates the validity of the basic approach. 

In conclusion, the impact of this inlet flow distortion analysis in the 
area of inlet/engine compatibility could be substantial. A means is now avail- 
able to carry-on compatibility studies prior to engine test. The analysis 
identifies those inlet flow and compressor design variables most important in 
the interaction of the inlet flow and the compressor. These are specifically 
the inlet distortion profile and magnitude, and the compressor reduced fre- 
quency. This latter parameter is directly dependent on the rotor chord which 
emerges as a strong factor in the design of hardware for compatibility. There- 
fore this approach can help to insure inlet/engine compatibility prior to hard- 
ware commitment. Furthermore, the successful treatment of the problem by 
fundamental aero-thermodynamic relationships can also put into perspective and 
tie together the number of empirical distortion indices currently in use. 

The favorable results of the data/analysis comparison are considered 
verification of the fundamental hypothesis of the analysis; specifically that 
distortion must be analyzed from an unsteady point of view. The accuracy is 
also considered sufficient to justify the simplifying assumptions incorporated 
in the analysis including the use of the overall compressor work balance rather 
than a detailed stage-by-stage analysis. Additional comparisons with a few test 
programs are recommended for refining the steady state analysis. The programs 
should provide information such as the detailed distortion patterns, compressor 
geometry and compressor operating conditions. However, the compressor model 
has been developed to the point necessary for conducting the phase of study 
designed to analyze the compressor reaction to unsteady turbulent flow. 



TASK II 

FLUID DYNAMIC MODEL OF TURBULENT INLET FLOW 


Turbulent flow produced in an aircraft inlet system can result in 
momentary total pressure distortion levels of a magnitude and duration 
sufficient to cause engine compressor stall. The first attempt to identify 
such instantaneous distortion levels was in terms of the Root Mean Square (RMS) 
level and the Power Spectral Density (PSD) function of the total pressure 
fluctuations. These statistical averages are relatively inexpensive to obtain. 
However, no physical interpretation could be given such quantities and as a 
result were of little value in the correlation of unsteady flow and compressor 
stall. Presently, it is common practice to measure an instantaneous distortion 
pattern at each instant of time by use of high response, highly (time) corre- 
lated total pressure instrumentation. This requires complex and expensive 
data measurement systems. While these instantaneous patterns graphically 
demonstrate the existence of unsteady flow, only limited empirical correlations 
of unsteady flow and compressor stall have been shown. The high cost of data 
reduction severely restricts the quantity of data analysis that can be made to 
substantiate any correlation. In addition, an empirical approach is inherently 
weak since it does not provide physical understanding of the basic flow phe- 
nomena. An analysis relating turbulent flow phenomena in the inlet to com- 
pressor stall is required. 

A fluid dynamic model of turbulent flow is developed herein as a means 
of understanding turbulent inlet flow and as a practical tool for evaluating 
flow properties through the use of the total pressure RMS level and PSD func- 
tions. By use of this flow model, which is based on a combination of basic 
fluid dynamic concepts and statistical analyses, a better understanding of the 
mechanics of the flow is obtained. Consequently, the loss in compressor stall 
margin may ultimately be related to the statistical characteristics of turbu- 
lence. 

The approach used to analyze this turbulent flow is outlined below and 
presented in detail in the subsequent sections of the report. Turbulence, 
normally measured in terms of velocity or total pressure, implies pressure 
gradients exist in the flow. It is a fundamental of fluid mechanics that 
pressure gradients can be supported by only two means: (l) pressure waves 

traveling at (or above) the local sound speed and (2) by streamline curvature. 
Since turbulence is produced by viscous phenomena, it is proposed that the 
pressure fluctuations measured in an inlet are primarily supported by stream- 
line curvature. To realistically model the flow, it is hypothesized that the 
streamline curvature and resultant pressure fluctuations are caused by a 
random distribution of discrete vortices being convected downstream by the 
mean flow. This is illustrated schematically in Figure 23. 

To obtain a mathematical representation of this flow model, a coordinate 
system moving downstream with the mean flow is used and enables the individual 
vortices to be analyzed in a steady frame of reference. Steady state flow 
equations can then be applied to describe the flow field about an isolated 
vortex. Analytical construction of the total pressure signature of this 
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isolated vortex is accomplished through the superposition of the vortex flow 
and the transport velocity. Consistent with the basic hypothesis, the turbu- 
lent nature of the flow is assumed to result from a distribution of these 
vortices having random size, strength, location, and direction of rotation. 

The total pressure root mean square level and power spectral density function 
resulting from this random flow field are found by use of statistical methods 
as applied to the analysis of this stochastic process. This development is 
followed by sensitivity studies to determine the impact of certain assumptions 
on the results. In addition, results' of the analysis are compared with test 
data. 



Pigure 23. Hypothealaed Turbulent Plow Composed of Random Vortices. 


Isolated Vortex 

The definition of the flow field associated with an isolated vortex is the 
first step in the development of this turbulent flow model. The description is 
based on a time dependent solution of the Navier-Stokes Equations. This vortex 
and associated properties are used in the subsequent development of the fluid 
flow model. 

Solutions of the Navier-Stokes Equations of Motion. - Depending upon the 
initial and boundary conditions, several vortex flow fields are found to be 
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solutions of the Navier-Stokes Equations. The details and the characteristics 
of these solutions are presented in Appendix F. It was found that all steady 
state solutions have singularities (infinite velocities) at either the vortex 
center or outer extremity and therefore do not represent real flows. There- 
fore, the time dependent solutions, each representing different vortex boundary 
conditions, were examined. For a realistic flow, the solution must satisfy the 
following boundary conditions: (l) the vortex must have a tangential velocity 

of zero at both the center and at an infinite radius, (2) the velocity must be 
continuous for all radii, and (3) the influence of the vortex must move outward 
with time. The solution of the two-dimensional time dependent Navier-Stokes 
Equations that fits these assumed boundary conditions is a vortex formed by an 
impulsive start in an undisturbed flow. The normalized velocity field asso- 
ciated with this vortex is given in Equation 22. 


■0 


©max 




( 22 ) 


where v^ 

v 0max 

a 


velocity in the 0 direction (cylindrical coordinate 
system) 

maximum velocity of the given vortex at a given time 
radius at which v Q = v^ 


n = specifies direction of rotation, = 1 for counter- 
clockwise and -1 for clockwise 


The radius r = a at which the velocity is a maximum is considered the vortex 
core radius. This core radius varies with time and as a result, the influence 
of this vortex increases in the radial direction with time. 


A complete development of this vortex is presented in Appendix H and 
includes a description of the angular momentum, vorticity, circulation and a 
discussion on the time of origin and decay of the vortex. 

Vortex Description in Cartesian Coordinates. - In order to describe the 
properties of the vortex in terms of a coordinate system fixed with respect to 
the inlet it is first necessary to express the vortex properties in terms of 
a Cartesian coordinate system fixed to the center of the vortex (x, y) as 
opposed to the cylindrical coordinate system (r, 0) used for solutions of the 
Navier-Stokes Equations. A description of the coordinate system used is shown 
in Figure 24(a). 

As developed in Appendix H the circumferential velocity is a function of 
the radius and time and is given by Equation 23. 

r 2 

Br ~ 4vt ( 2 3) 

ve - pr e 
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■where: B = constant depending on the vortex strength 


l> = kinematic viscosity 
t = time since the vortex started. 


However, due to the short period of time the vortex is in the field of interest 
(the inlet duct) and the very slow rate of growth of the vortex, the time is 
assumed constant. The vortex tangential velocity normalized by the maximum 
velocity, which occurs at r = a, was given by Equation 22 and repeated below. 
The variation of Vq with r is shown in Figure 24(b). 


-^[(r/a) 2 -!] 

V0 - n(r/a) e 

v Omax 


( 22 ) 


The horizontal (u) and vertical (v) velocity components of the vortex 
velocity are obtained by use of the description and definitions of the coordi- 
nate systems of Figure 24(a). These are: 


->5 [(x/a) 2 +(y/a) 2 -l] 

u ” - v 6ma5 (Y/a) e 


(24) 


[(x/a) 2 +(y/a) 2 -l ] 

v " v 0mai< x/a) 6 


(25) 


The local flow angle is dependent on the velocity components, and is given by 
Equation 26. 

a ■ arctan (v/u) ■ arctan ( tan ^g) (26) 

The static pressure variation with radius as determined in Appendix I is given 
by Equation 27 and shown in Figure 24(c). 

« „ _-q „2 [ (x/a) 2 +(y/ a) 2 -l] 

static pressure at radius r 

static pressure beyond the vortex influence, 
r > > a 

density 

The total pressure in the vortex flow field, based on a coordinate system fixed 
to the vortex center, is equal to the sum of static and dynamic pressures or 


r r ° " 2 ®max 


where; 


P o = 


P = 


4i 



(28) 


P T 


Pr+f v 6 


2 


In terms of the cartesian coordinate system, the total pressure is obtained by 
combining Equations 22, 27, and 28 and is given by Equation 29. 


P T -P o " 2 [ (x/a)2+( y /a)2 - 1 ] v 9max 6 


[(x/a) 2 +(y/a) 2 -l] 


(29) 


Transformation of the Vortex Flow Field to the Inlet Coordinate System. - 
The vortex flow field has been described in the coordinate system fixed to the 
vortex center. The vortex, however, is moving downstream at the local average 
flow velocity. To determine the total pressure as it would be measured at the 
engine face, it is necessary to transform the vortex velocity field into the 
coordinate system fixed to the inlet. 

Allowing the vortex to move downstream at the local velocity, U Q , as shown 
in Figure 25, the local instantaneous velocity components as would be measured 
in a coordinate system fixed at the probe becomes: 


U “ V^^/a) e 


$i[(X/a) 2 +(Y/a) 2 -l] 


(30) 


V “ ” V ®max (X/a) 6 


h [ (X/ a) 2 +(Y/ a) 2 -l] 


(31) 


Here, the upper case (X, Y) represent the position of the vortex center relative 
to the origin of the fixed coordinate system located on the probe. In essence, 

U and V are the velocity components at the probe due to a vortex located at 
(X, Y). In this case, the transformation between the moving and fixed coordi- 
nate systems is simply 

X = 

Y = 

The static pressure in the fixed coordinate system is 

P „. Po . _£L ,7 e - I(X/a) 2 +(V/a)2-l] (33) 

r r o 2 e max 


(32) 

y 
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Figure 25 • Transformation of Coordinate Systems. 


To account for the fact that the local flow angle at the total pressure probe 
is not aligned with the probe, the total pressure as expressed in Equation 34 
is assumed to be the local static pressure plus a corrected dynamic pressure. 


PT - P r +n R (a)|w 2 


(34) 


where: P = static pressure (which is independent of coordinate 

r system) 

W = the resultant velocity vector, (U 2 + v 2 ) 1 / 2 
tj = probe total pressure recovery 

XV 

a = local flow angle 

Since the vortex is moving downstream at velocity IL, the position, X, of 
the vortex center with respect to the fixed (probe) coordinate system is a 
function of the velocity, U 0 , and time, t. Thus: 


X - u 0 t 


(35) 


This implies that the origin of time must be chosen such that X = 0 when t = 0. 

The total pressure recovery is assumed to vary directly by the square of 
the cosine of the flow angle, a , at the probe. 


2 ,U N 2 

n R “ co8 z a - (-) 


(36) 
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Using the transformation of Equation 35 , the total pressure recovery, and 
the definition of static pressure (Equation 33 ) and velocity (Equation 30 and 
31), the total pressure expressed by Equation 34 becomes: 


P T -P T - n pU 0 vq 


max a 


# 2 -u 

(-)e 


( 37 ) 


+n 




max 


[(V/a) 2 -l 


_[(Uat) 2 +( X) 2 _i] 


where P To 
P 0 + o 


is the total pressure of the undisturbed flow and equal to 
2 


Defining this delta pressure as A P T and normalizing by the dynamic pressure, 


<Jo 


* U„ 


( 38 ) 


the total pressure becomes: 

ap v 0_ -!:[(U 0 t/a) 2 +(Y/a) 2 -l] 

APx . 2 n(-^2L) (Y/a) e 

u o 


+n 2 (liBSS) 
V U rt ' 


[(Y/a) 2 - 1 ] 


- [(U 0 t/a) 2 +(Y/a) 2 -l] 


( 39 ) 


The resultant total pressure is a function of the time, t, the mean 
velocity, U , the strength of the vortex VQ max , the core size of the vortex 
a, and direction of rotation, n. Detailed development of these relation- 
ships are given in Appendix I. 

Although these relations have been developed for a single vortex, they are 
applicable to a series of vortices. These vortices can assume random values 
of size, strength, location and direction. In such a case the vortex proper- 
ties, a, y, and n, become random variables. The techniques required 

to treat these properties as random variables are developed in Appendix J . 

These tools are then used in conjunction with Equation 39 to establish the 
statistical characteristics (RMS and PSD) of the proposed flow field. 


Statistical Flow Model 

The model of turbulent inlet flow is hypothesized as being composed of 
a random distribution of vortices, each having a specific size, strength, 
direction of rotation, and location as shown in Figure 23 . The total pressure 
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fluctuation created by each vortex is given in Equation 39* For a specific 
vortex having a given set of properties, a, v p^y ? Y> n > Equation 39 
signifies a single time function. However, each vortex has a different set of 
properties; consequently, the flow field is composed of a family of time 
functions. This family is called a stochastic process. 

The autocorrelation function and its Fourier Transform, the power spectral 
density function resulting from this process, are found in functional form by 
statistical methods as applied to a general stochastic process in Appendix J . 
These developments are applied to the vortex flow field to obtain the total 
pressure autocorrelation and power spectral density functions of the turbulence . 
The fluctuations in total pressure are also related to the fluctuations in 
velocity by use of the flow model, providing an analytical means to relate, 
hot wire anemometer velocity measurements and total pressure measurements. 

Autocorrelation Function. - The autocorrelation function of the stochastic 
process composed of the random vortices flowing downstream with the flux of N 
per second cam be summarized as follows. 

The general total pressure wave described by Equation 39 can be repre- 
sented by the following functional notation. 

AP t = AP t (a,v,Y,n,t) (40) 

The autocorrelation function of this discrete total pressure wave is found by 
means of the definition of the autocorrelation function and is given by 

oo 

R (a,v, Y,n,t) = / AP T (a,v,Y,n,t) AP ? (a,v,Y,n,t + t)dt (4l) 

APt J l 


To establish the autocorrelation function of the entire family of waves, a 
weighted sum of the function (a, v, Y, n, r) is performed over all 

possible values of a, v, Y, and n. T The weighting functions are the percentage 
of vortices having a specific core size between a and a +Aa, specific strength 
between vq,^ and v^^ + Av^, a specific location between y and y +Ay and 
a direction of rotation (n = +1 or -l). These weighting functions are assumed 
independent of' each other, are simply the probability density functions of a, v, 
Y, and n, respectively, and are designated by the notation, p( ). The auto- 
correlation function of the resultant total pressure wave (Equation 42) is there- 
fore the weighted sum of the general autocorrelation function. 


V <T) 


N *W T a » v » Y » n » ' T ) p ( ▼) ' p ( Y) P( n)dndYdvda 


(42) 


The autocorrelation function of the vortex flow field, as measured at the 
total pressure probe, is found by incorporating the total pressure wave 
(Equation 40), the definition of the autocorrelation for a discrete wave 





(Equation 4l), and appropriate probability density functions into Equation 42. 
This substitution is accomplished in Appendix J and the resulting integral 
equation is solved by numerical techniques. The particular probability density 
functions used in this analysis are discussed in detail in Appendix K and out- 
lined in Figure 26. 

A computer program was written to evaluate the integral and the resulting 
normalized autocorrelation function is shown in Figure 27. Although an abso- 
lute autocorrelation function is directly dependent on the probability density 
functions of the mean core size, a, and the mean vortex strength, v^g^, it 
will be shown in the sensitivity studies this normalized function shown in 
Figure 27 is relatively insensitive, to the probability density functions used. 

RMS Total Pressure Fluctuations: The mean square total pressure fluctu- 

ation is identical to the autocorrelation function at a time delay, T , equal 
to zero. Thus, 

(7 2 ■ APp ■ !?AP T (o) (43) 


For this case and using the assumed probability density functions (Figure 26) 
with h = (1/2 )H (probe at center of the duct) and a/H«l Equation 42 can be 
integrated in closed form to give 



lle^ (jBjjr+2_)__ ( m v + 3 ) (m v +U ) 

32 (m v +n v +3) (m v +n v +4) (my+ny+5) . 


The mean square level is later used to normalize the power spectral density 
function and utilized below to relate the velocity and total pressure fluctu- 
ations. 

Relationship Between Velocity and Pressure Fluctuations: The flow model 

also provides a base from which other quantitative relationships can be 
developed. By application of the statistical techniques developed herein, the 
fluctuating velocity components defined by the isolatedvprtex, the auto and 
cross correlations of the unsteady velocity terms (u^, v^, uv) can be obtained. 
The detailed development of these correlations are presented in Appendix L. 
Briefly, the mean square axial perturbation velocity in the center of the 
duct is: 

ii 2 eim 2 ,NH* .a. , V 9max > (m a +2) (m^+2) 

SP “ 2 MJ 0 <..V ' U c ' (m a +n a +3) (*^+3) (45) 
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The relationship between the velocity fluctuations and the total pressure 
fluctuations in a turbulent flow field is obtained by combination of Equations 
44 and 45. This is: 


U 2 /Ur > 2 L_ 

APT Z /qo 2 . . lie (my+3) («v+4)_ 

' 16 (tftv+nv+4) (my+ny+5 


(46) 


For the velocity probability density function having the exponents niy = 4 and 
n v = 14 (Refer to Figure 26 for the density functions) the ratio has a value of 


= 0.238 (47) 

^Wo 

It should be noted that for reasonable values of % and n^., the numerical value 
of Equation 46 is weakly dependent upon these variables. 

Equation 47 can be written in terms of the mean Mach number, M 0 , and total 
pressure, Pip, as follows and is shown graphically in Figure 28. 

~2 Ap 2 2(1+^M 2 ) y_1 

« 0.238 — ^ 

This result, in itself, is significant. For the first time a relationship has 
been developed which will relate the fluctuating velocity as measured by hot 
wire an emom etry with the total pressure fluctuations. Previously, the relation 
ship was developed by ass umi ng either sonic waves or a quasi-steady analysis 
with a constant static pressure, neither of which represented the physical 
process. 

Power Spectral Density Function. - The power spectral density function is 
the Fourier Transform of the autocorrelation. Thus: 



S (f) - (T)e" 32 ' £T dT 


■viiere: S AP T ^ f ^ = 

R AP t (t) = the 


(48) 

complex power spectral density function, 
autocorrelation function. 
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Figure 2& Relation Between Velocity and Pressure Turbulence. 

Of immediate interest is the real part which is the physically realizable 
one-sided power spectral density function as normally obtained from test 
data. This function G (f) is defined (see, for example, Reference 11) as: 

/*00 

<U f > " 4 / R Mv (T)cos(2lTfT)dT (49) 

0 

The normalized power spectral density function obtained from the auto- 
correlation function illustrated in Figure 27 is shown in Figure 29. The 
spectrum is normalized by the mean square total pressure fluctuations, c r^, and 
the mean vortex core size, a. 

Sensitivity Studies 

The effect of the assumed statistical distributions on the resultant 
autocorrelation and power spectral density functions must be determined since 
these distributions, in effect, describe the mean core size and mean strength. 
The level and frequency content of the power spectral density function is 
directly related to the mean vortex strength as measured by vg^ax/lJ * the mean 
core size, a/H, and the shape of their respective probability density curves. 
However, the normalized spectrum may not be sensitive to these variables, since 
it is normalized by the mean square of the pressure fluctuations and the mean 
core size. If this is the case, a great simplification in application will 
result. It is the specific objective of this section to establish this sensi- 
tivity. 

While it was felt that the power spectral density function is primarily 
dependent on the strength and size of the vortices, the effects of spin direc- 
tion, n, and probe location, Y/H, were also investigated. The effects of these 
variables on the power spectrum were found to be negligible. 
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NOTE: The probability density functions 

are given in Figure 26. 



In the turbulent flow model, the vortex strength is represented by the 
maximum tangential velocity. The distribution of vortices having various 
strengths is governed by a Beta probability density function given in Equation 
50: 


p. .. k v / v 9 m ax V 
V ®max " U Q ( U Q ) 


\ m v 


1 - 


V 9max\ nv 

U„ 


(50) 


where 


max 


U 


m +1 
v 


m + n +2 
v v 


and 


K 


T (m + n + 2 ) 
v v 

r(m + 1) r(n + 1) 
v v 


U V V V V 

The shape of this density function can be altered by changing the constants 
ny and n v as illustrated in Figure 30. The effect of such changes serve only 
to modify the mean square value of the fluctuating pressures and do not affect 
the shape of the normalized autocorrelation or power spectral density function 
as shown in Figure 31. 



Figure 30* Probability Density Function of Vortex Strength. 

The distribution of vortices having various core sizes is also described 
by a Beta probability density function as given by Equation 51: 


P( a) 


ka /aA m a 
H H 


where 


a. 

H 


1-5 


m +1 
a 

m + n + 2 
a a 


K 


r(m + n + 2) 
a a 

a = f(m a + 1) f(n a + 1) 


( 51 ) 
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NORMALIZED AUTOCORRELATION - R APi (t) / (? aPt (0) 



NORMALIZED TIME DELAY - TUo / a 


Figure 31 (a). The Effect of Velocity Density 
Function on the Autocovariance 
Function 


NORMALIZED POWER SPECTRAL DENSITY 



0.01 0.1 1.0 

NORMALIZED FREQUENCY - f a/EJ 0 


Figure 31(b). The Effect of Vortex Velocity Density Function on 
the Normalized Power Spectrum. 


Variations in shape of this function have some impact on both the amplitude of 
turbulence and the autocorrelations and power spectral density functions. 
Several variations used in the sensitivity study are shown in Figure 32. An 
example of the effect of these variations on the autocorrelation and power 
spectral density function is shown in Figure 33. 



Pigure 30 ..Probability Density Function of Core Size. 


As illustrated in Figures 3l("b ) and 33(b), the normalized power spectral* density 
function is independent of the mean value of the vortex strength v @ m ax and the 
shape of its probability density function and only weakly dependent upon the 
mean vortex core size a and its probability density function. If these latter 
variations are ignored for the moment, it is apparent that this normalized 
spectrum along with particular values of the mean core size a and the mean 
square of the pressure fluctuations cr^ can ^e used to calculate an absolute 
PSD function (G^p^C^) as a function of f) .Furthermore, if the power spectrum 

is normalized by cr^, the resulting absolute spectrum is only a function of a. 

The importance of this result lies in the fact that the analytical spectrum 
can be easily matched to experimental data using only the normalized spectrum 
and the mean core size a. As an example of this concept the analytical model 
was matched to the turbulence data from an axisymmetric mixed compression inlet 
presented in Reference 12. The procedure for determining the mean core size 
associated with this particular spectrum is given below and illustrated in 
Figure 3^. 

1) Normalize the measured power spectral density function by the mean 
square of the turbulence ( 0 ^). 

2) Compute an absolute spectra (G from the normalized 

spectrum obtained from the analysis by assuming various values of the 
mean core size, a, and using the local velocity, U Q . 
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NORMALIZED AUTOCORRELATION - R ap /t) / R AP ^0) 



NORMALIZED TIME DELAY - TUo / a 

Figure 33(a). The Effect of Vortex Core Size 
Density Function on the Auto- 
covariance Function, 


NORMALIZED POWER SPECTRAL DENSITY 


NORMALIZED FREQUENCY - f a/U 0 

Figure 33(h), The Effect of Vortex Core Size Density Function on 
the Normalized Power Spectrum. 



3) Compare the computed spectra with that obtained experimentally. The 
mean core size producing the best fit with the data corresponds the 
characteristic mean core size for the particular flow conditions. 

As shown in Figure 34, the mean vortex core size producing the best fit under 
the constraints of both the shape of the spectrum and the area under the spec- 
trum (cr 2 ) is 2.1 centimeters. The comparison shown yields excellent agree- 
ment. Additional examples of test/analysis comparison will be discussed in the 
"Data/Analysis Comparison" section. 

As previously indicated and illustrated in Figure 33(b) the mean core size 
and shane of its orobability density function have a small effect on the normal- 
ized power spectral density function. Therefore, depending upon. the particular 
spectrum used, this frequency shift in the normalized spectrum will result in. 
variations in the mean core size, a, obtained from the analytical model when it 
is matched to experimental data. As a measure of the sensitivity of the com- 
puted mean core size to the assumed probability density function, the computa- 
tions summarized in Table V were made. It was assumed that the test data was 
represented by the normalized power spectrum with the core size and maximum 
velocity probability density functions corresponding to = 2, n a = 2b, and 
™ = 4, n = 14, respectively. This result is indicated by the asterisk. m 
Table V and shown in Figure 32. The frequency shift of the other normalized 
spectra, caused by other assumed density functions, will thus cause a false 
indication in actual mean core size. This error is given in Table V as A a/a*. 
As shown the maximum error is 25$ and occurs for assumed probability density 
functions yielding large values of a/H. For data indicating mean core sizes 
greater than 30$ it may be necessary to choose coefficients that give consist- 
ent values of core size. This would require iteration. However, it is felt 
that tynical values of core size will be less than 30$, in which case the error 
will be" small and hence iterations will not be necessary. 

Results of this sensitivity study indicate that for turbulent flow with 
a./H < 0.30 the probability density functions of core size and maximum velocity 
have only small impact on the overall turbulent flow model itself. This is. 
important to the use of the model eliminating two degrees of freedom that might 
ordinarily have to be taken into account. 
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POWER SPECTRAL DENSITY 



FREQUENCY - f, H* 

Figure 34. Comparison of Spectra Computed from the Flow Model with 
Inlet Test Data of Martin (Reference 12) . 
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TABLE V 

IMPACT OF ASSUMED DISTRIBUTIONS ON POWER SPECTRUM 


m a 

n a 

a/H 

m v 

n v 

^6 msux/^o 

A a/a* 

2 

26 

.100 

2 

8 

.250 

0 

2 

26 

.100 

6 

20 

.250 

0 

2 

26 

.100 

8 

26 

.250 

0 

0 

8 

.100 

4 

14 

.250 

-.14 

1 

17 

.100 

4 

14 

.250 

-.06 

*2 

26 

.100 

4 

14 

.250 

0 

3 

35 

.100 

4 

14 

.250 

.025 

4 

44 

.100 

4 

14 

.250 

.055 

2 

2 

.500 

4 

14 

.250 

.245 

6 

2 

.700 

4 

14 

.250 

.173 

6 

14 

.318 

4 

14 

.250 

.145 

2 

26 

.100 

k 

34 

.125 

-.003 

2 

26 

.100 

5 

9 

.375 

.021 

2 

26 

.100 

4 

4 

.500 

.056 

10 

10 

.500 

4 

14 

.250 

.146 

20 

20 

.500 

4 

14 

.250 

.145 
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Scaling Law For Turbulent Flow 

Inlet development universally begins with subscale testing. These test 
results are then extrapolated to full scale to establish the expected inlet 
performance. The technique to scale turbulence has, however, not yet been 
defined, although dimensional analyses suggest that the power spectrum will 
shift in frequency in inverse proportion to the actual scale size. 

The scaling techniques have not been defined because of the complex nature 
of random, unsteady flow. "Ever since the derivation of the momentum equation, 
the fundamental problem of the analysis of turbulent flow has been that of 
closing the system of governing equations. This is caused by the fact that 
even in its simplest form the turbulent flow momentum equation contains a 
"Reynolds stress" term made up of a correlation of the fluctuating components 
of the turbulent velocity field which acts as an apparent stress. Since the 
momentum equation is the governing equation for the mean velocity field, the 
presence of this apparent stress term introduces additional unknowns into the 
problem, and any equation derived without further assumptions to characterize 
these unknown quantities will in turn introduce other unknown quantities. One 
method of closing the system of equations is to formulate models for the tur- 
bulent shear stress in terms of already known (or knowable) quantities 
(Reference 20). This formulation of models has classically been handled in one 
of two ways — either some model can be postulated for the turbulent shear stress 
itself, or, in analogy to a laminar flow, the turbulent shear stress can be. 
assumed to be given by some effective viscosity multiplied by a local velocity 
gradient . 

An attempt to develop a third, perhaps more fundamental model, was made 
by use of the fluid dynamic model of turbulent flow developed in this program. 
The model was used to define the "Reynolds stresses" in terms of the mean 
turbulent flow properties, providing a third model for the shear stress. 


This development is outlined in Part A below. The resultant set of 
Equations could not, however, be solved for an arbitrary inlet flow profile 
within the scope of the planned effort. Nonetheless, a specialized case is 
given in Part B which may give some indication of the method of scaling the 
turbulence power spectrum. 


A. General Development: The longitudinal velocity is assumed to be 

dominant in the following development and the flow velocity and static pressure 
is defined in terms of a fixed value plus the perturbation value as indicated 
below: 


V = 0 + v 

P - P Q + A P 

The initial velocity distribution in the x-direction serves as a boundary 
condition. This distribution is shown sketched below: 


(52) 
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Utilizing the notation defined above the following ten Equations can be identi 
fied. 

Momentum: 


[ 

3u — _3u_ 

= -I 

9u 

8uv 

0 

9x + v 9y 

p 9x 

9x 

' ” 3y 

1 

ll + v 

_ -1 9A£ _ 

9uv 

97 

0 

9x 9y 

p ay 

3x 

3y 


(53) 




Continuity: 

3]L_ + = (55) 

9x . 9y 

Crocco's Theorem: As applied to a velocity change caused by an entropy gradient 

in total pressure, Crocco's Theorem (See for example Reference 21, page 2ol) 
relates the total vorticity, NT , to the velocity gradient AU . This is: 
is : 

m Udu - u£ 2 - 0^ a B 4 B (56) 

2 

Fluid Flow Model: The remaining six equations are defined by use of the fluid 

flow model developed in this program. The alternative Equation numbers refer 
to their explicit definition as established elsewhere in this report. 
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( 57 ) 


u 2 = f (N, a, v 0 ) 

max 


V 2 = f(N, a, v 0 ) ( 58 ) 

max 

uv = f(N, a, V- ) (59) 

0 

max 

u = f (n, N, 7 , Vg ) (6°) 

max 

v = f(n, N, a, v 0 ) (^l) 

max 

AP = f(N, a, v ) (62) 

0 

max 


Similarly, the following ten unknowns are identified: 


u, v, Ap, u , 


2 — 
v , uv. 


N, 


n 


a, 



max 


(63) 


By solution of these ten equations, the unknowns u^ and a can be obtained as a 
function of the initial velocity (distortion) profile. These quantities are 
directly related to the total pressure RMS level and power spectrum as estab- 
lished by Equation 46 and the techniques leading to the spectrum of Figure 29. 
If it can be assumed that the inlet velocity (or total pressure) profile can be 
scaled, the techniques for scaling RMS level and power spectral density func- 
tion will result. 


Solution of these equations for an arbitrary distortion profile was beyond 
the scope of the planned program. However, some insight can be gained into 
scaling of the power spectrum by application' to a specialized case. 

B. Specialized Velocity Profile: Consider flow in a two dimensional duct 

that has a non-uniform total pressure and hence velocity profile. As estab- 
lished from fluid flow model, a gradient in velocity will occur if the vortices 
have a non-uniform distribution of direction of rotation, e.g., more vortices 
having a positive spin than negative spin. The general expression for this 
gradient is given by Equation L-13 in Appendix L. 
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Equation L-13 can be integrated in closed form at the duct walls, h = 0 and 
h = H, by assuming a/H«l. This results in 


u 

U„ 


h=H 


0 — 2 

,NH n s 1/2 , max N ,a s 

= -n (— ) ( 27 re) (— ) (-) 

U 0 U 0 H 


(64) 
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(65) 


This will enable the maximum difference in velocity across the duct to 
be established. The total vorticity flux is dependent only on this differ- 
ence and hence for the purposes of this example it is not nexessary to obtain 
the velocity profile between h = 0 and H, which would require numerical inte- 
gration of Equation L-13. However, it is for this special case of interest 
that the scaling technique will be defined which is valid only for the velocity 
profile implied in Equation L-13. 


The absolute value of vorticity passing between h = 0 and H can be com- 
puted by use of the flow model as in Appendix H. This flux of vorticity 
becomes : 

— v e 

Nr = n 8tt NH U e" 1/2 (§) (-~ 2L ) (66) 

U H U 0 

However, this flux is related to the change in velocity across the duct by 
Crocco's Theorem, Equation 56, which in turn is obtained by combining Equations 
64, 65, and 66. By so doing it is found that a/H is a constant (within the 
confines of the velocity profile described by Equation L-13). 

Since the power spectral density function is directly related to the 
mean vortex core size, this result specifies that increasing inlet size would 
shift the spectrum in inverse proportion to the inlet scale. This result would 
verify the present practice of scaling the frequency content of the turbulence 
power spectrum inversely with the inlet scale size. Again, this applies only 
to the specialized case. 
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Data/Analysis Comparison 


The power spectral density functions of turbulent flow computed from this 
turbulent flow model were compared with data obtained from Reference 12 for a 
test of a .55 meter diameter, Mach 3 mixed compression axisymmetric inlet. The 
comparisons were performed to illustrate the technique of determining the mean 
turbulence scale size, verify the assumptions incorporated in the turbulent 
flow model, and to gain insight into production of turbulence by interpretation 
of test data through use of the model. 

The experimentally determined power spectra used in the comparisons are 
shown in Figures 35 and 38. The mean core size of the turbulent flow was com- 
puted for these spectra by assuming a mean core size and comparing the test 
data with results of the turbulent flow model. These comparisons are shown 
in Figures 36 and 37 and Figures 39 through 44 . Excellent agreement is 
obtained indicating the turbulent flow model accurately represents the real 
case and that the flow can be modeled by the assumption of random vortices with 
the indicated mean core size. 

Data from probe number 900 (Figure 36), which is in the inlet duct exten- 
sion, approximately 3-1/2 diameters downstream of the inlet diffuser exit 
(Figure 36) shows some disagreement with the analysis at the higher frequencies. 
However, if it is assumed that this scatter is created by turbulence from two 
separate sources each having a different mean core size, a composite spectrum 
can be constructed incorporating the two different mean core sizes. This com- 
posite spectrum is shown compared with data in Figure 45. Excellent agreement 
is obtained. This not only adds credence to the turbulent flow model but allows 
concrete interpretation of the power spectra. For example, it can be postulated 
that the source of the turbulence having the scale size of approximately 
a = 2.0 cm. at the compressor face probe 872 may be produced by shock boundary 
layer interaction. This turbulence carries through to the inlet duct extension 
(probe 900) but is reduced in size to approximately a = 0.8 cm., as indicated 
in the composite of Figure 45. However, a new source is identified by intro- 
ducing of turbulence having the scale size of 2.9 cm. which may be attributed 
to the wakes generated by the centerbody support struts. 

As a second, perhaps more significant example, consider the data of Figure 
39 (Probe 872) in which case the total pressure recovery was 88$ and the mean 
core size of the turbulent flow 2.14 cm. As the inlet is operated at increased 
supercritical margins (increased terminal shock strength) the power spectrum 
changes as shown by the data in Figure 44 for a recovery of 57$. 'A composite 
spectrum for this supercritical case was computed from the turbulent flow model 
and is shown in Figure 46. It indicates the flow is composed of eddies having 
a mean core size of 1.9 cm plus additional turbulence having a mean core size 
of 5.0 cm. This may be indicating large scale separation is occurring and the 
turbulence is not only increasing in strength but in size as well. 

Comparison between the power spectrum as established from the analytical 
model and that obtained from test data shows exceptional agreement giving 
considerable credence to the development outlined herein. With this model the 
total pressure power spectrum and root mean square level of the total pressure 
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fluctuations take on considerable significance. It is specifically the 
strength and size of these low pressure regions, derived by application of the 
model, that are important in the inlet flow/ engine interactions. 
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Figure 35 - Comparison of engine face and duct exit total pressure power spectral densities, 
Mq = 3.0. (Reference 12) 
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Figure 37. Comparison of Spectra Computed from the Flow Model with 
Inlet Test Data of Martin (Reference 12). 
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Figure 39. Comparison of Spectra Computed from the Flow Model with 
Inlet Test Data of Martin (Reference 12) . 
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Figure 40. Comparison of Spectra Computed from the Flow Model with 
Inlet Test Data of Martin (Reference 12) . 
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Comparison of Spectra Computed from the Flow Model with 
Inlet Test Data of Martin (Reference 12) . 
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Figure 43. Comparison of Spectra Computed from the Flow Model with 
Inlet Test Data of Martin (Reference 12). 
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Figure 45. Construction of the Composite Spectrum and Comparison 
with Data from Probe 900 (Downstream) 
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CONCLUSIONS AND RECOMMENDATIONS 


The an al ytical developments reported herein provide a fundamental approach 
to the problem of inlet/engine compatibility. With further development, this 
approach will provide a method of evaluating inlet tests and engine designs 
early in the propulsion system development cycle. Ultimately, it shows promise 
as a method for predicting and evaluating the effects of distortion and turbu- 
lence on engine stall characteristics, prior to system test. 

The following are the more significant conclusions arising from the work 
to date and suggestions for continued activity to achieve the basic program 
goal: establishing the fundamental relationship between inlet distortion and 

turbulence and loss in compressor stall margin. 

Conclusions 

(1) The effect of circumferential total pressure distortion on the loss 
in compressor stall margin has been established analytically. The analysis has 
shown that the stall margin loss_is directly a function of the distortion 
pattern, the distortion level ((P t - Pt min )A)> and of the compressor rotor 
reduced frequency, k =u c/ 2 U. 

(2) The rotor chord is the principal design variable in the reduced 
frequency, k, and therefore emerges as a significant engine parameter in 
design for compatibility. 

( 3 ) Favorable comparison of distortion and engine stall data with analysis 
results is considered verification of the fundamental hypothesis of the 
analysis. Specifically, a circumferential total pressure distortion will result 
in an unsteady flow over the rotor blades requiring these unsteady aerodynamic 
effects to be included in the stage characteristics. 

(4) The accuracy of the stall prediction technique is sufficient to 
justify the simplified approach which considers an overall compressor work 
balance rather than a detailed stage-by-stage development. 

( 5 ) A phenomenological model of turbulent flow typical of that found in 
aircraft inlets has been developed by combining statistical techniques with 
the basic laws governing fluid dynamics. The power spectral density function 
and root mean square level of the fluctuating total pressure take on consider- 
able significance as a consequence of the model resulting in a means, of 
determining the strength and extent of time variant low pressure regions. 

( 6 ) Favorable comparison of spectra obtained from the analytical model 
with test data of a Mach 3 mixed compression inlet verify the Turbulent Flow 
Model. 

( 7 ) The ag reement with test data for both the Compressor Analysis and 
Turbulent Flow Model strongly suggest that compatibility problems, heretofore 
only attacked by empirical methods, are amenable to analysis. 
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Rec ommendations 


(1) Both the Compressor Analysis and Turbulent Flow Model are considered 
developed to the point necessary to initiate the program to achieve the long 
term goal of establishing a fundamental relationship between both inlet dis- 
tortion and turbulence and compressor stall margin loss. 

(2) Further comparisons of the compressor analysis with test data are 
recommended for refining the method. The data used should provide the detailed 
distortion patterns, compressor geometry and compressor operating conditions. 

(3) Additional analysis of turbulence data from a well documented inlet 
test program should be conducted to demonstrate the use of the turbulent flow 
model in isolating the source of turbulence and establishing turbulence decay 
characteristics . 

( 4 ) Finally, the developed relationships between distortion, turbulence 
and loss in compressor stall margin should be compared with data from an 
aircraft flight test program to verify the analysis. 
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APPENDIX A 


Analysis of Unsteady Potential Flow on an Airfoil 


This appendix contains a detailed discussion of the physical mechanisms 
effecting the interaction between an airfoil and the surrounding unsteady flow. 
The effects of unsteady airflow or airflow motion on an airfoil below the 
stalling angle of attack are explained in terms of the physical mechanism in 
part A and then related quantitatively with the governing mathematical rela- 
tionships in part B. 

A. The Governing Physical Mechanisms . - Lift on an airfoil is a conse- 
quence of the unequal pressures acting on the upper and lower surfaces. In 
potential flow these pressures along the surface of an airfoil can be computed 
from the velocity field by use of the equations of motion. The result of a 
potential flow solution, with zero circulation, about a flat plate of length, 
"c", and inclined to the stream at an angle, a , would yield stream lines similar 
to those sketched in Figure Al. 



Figure Al: Potential Streamlines Around an 

Airfoil with no Circulation 


The required boundary condition is that of zero flow through the plate. This 
is equivalent to the plate being a streamline which by definition implies zero 
cross flow. The lift on the airfoil is given by Equation A-l and is designated 
as the n on-circulatory lift, Although the lift component is zero for 

steady flow, it can be non-zero under non-steady conditions. 

l NC m ~ o < P U ” p l) dx ( A_1 ) 

To bring this mathematical model into agreement with experience it is 
necessary to move the rear stagnation point aft to the airfoil trailing edge 
(Kutta condition). This can be accomplished by imposing a circulatory flow 
around the airfoil as shown in Figure A2. 
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NON-CIRCULATORY FLOW PURE CIRCULATORY FLOW COMBINED FLOW 

Figure A-2. Flow About an Airfoil with Circulation. 


By so doing an equal circulation of opposite sign is shed from the airfoil such 
that the sum of the two is zero. The circulation about the airfoil imparts an 
additional velocity field around the airfoil and must be accounted for in 
obtaining the total lift. This additional lift is directly proportional to the 
amount of circulation r required to move the stagnation point aft to the 
trailing edge and is given by Equation A-2. 


l c - pur 

This in turn is related to the angle of attack, a , of the airfoil by 
“ pU o r - 2ira(-|pU 0 2 )c 


(A-2) 


(A-3) 


where r ■ irU 0 ca 


(A-4) 


Total lift is therefore composed of two terms each caused by imposing a boundary 
condition: 

L - L N c + L c (A-5) 

where is of non-circulatory (potential) origin and is a consequence 
of the requirement for zero flow through the plate 

and L c results from the circulation required (Kutta condition) to move 
the rear stagnation point to the trailing edge. 

In steady flow and as indicated by the symmetry of Figure A-l, the lift 
due to non-circulatory flow is zero. Thus, in steady flow the total lift is 
approximated by the circulatory term, 

L - L C - (^pU 2 )c 2im (A-6) 
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For unsteady airfoil motion, however, inertial forces due to the finite 
mass of fluid that is forced to adjust cause the non-circulatory lift to be 
non-zero. As an example, consider the resultant lift due to oscillatory 
vertical translation and pitching motion (about the mid chord) of an airfoil: 

L nc ■ irp ^ [y + Ua] (A-7) 

c 2 

where irp -jj- is termed the virtual mass per unit span associated with the 
vertical acceleration y and/or Ua. It is equivalent to the mass of a cylinder 
of fluid with a diameter equal to the airfoil chord. 

In addition, for any change in motion of the airfoil a trailing edge 
vortex must be shed to force the stagnation point to remain at the trailing 
edge. This adjustment requires a finite time and it too contributes to the 
resultant airfoil lift. The vortex shed is non-recoverable and might be 
considered as a dissipative or damping force. For the case mentioned above 
the lift of circulatory origin is: 

1>C ■ TrpUcC(k) [ y + Ua + ^ a ] (A-8) 

Combining equations A-7 and A-8 the total lift for an airfoil tindergoing 
oscillations in the vertical direction and along the pitch axis is as follows 

L - irp f 2 [ y + Ua ] + irpUcC (k) [ y + Ua + f «] (A-9) 

For a general treatment of non-steady airfoil/airflow refer to References 
18 and 19. 

B. The Governing Equations . - The pressure forces due to unsteady flow 
acting on the isolated airflow discussed in Part A will be quantitatively 
described below. This will be important in developing the effects of unsteady 
flow on the stalling airfoil lift coefficients in subsequent sections. 

Airfoil lift in unsteady airflow is dependent upon three prime factors: 

1. The type of motion occurring. This can be classified by the 
following : 

(a) Airfoil unsteady - freestream airflow steady 

(i) Airfoil undergoing vertical translation 

(ii) Airfoil undergoing horizontal translation 

(iii) Airfoil exhibiting pitch oscillation 

(b) Free stream airflow unsteady - airfoil steady 

(i) Airflow with vertical gust 
(ii) Airflow with horizontal gust 

2. The type of input, e.g., sinusoidal, step, arbitrary, etc. ‘ 

3. The times associated with the disturbance compared with the airfoil 
"time constant". 

Unsteady lift equations for the various classes of motions are summarized in 
Table A-l. 
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TABLE A1 UNSTEADY LIFT EQUATIONS 


MOTION 


SKETCH OF INPUT 


(0) Steady State 


% 

etc. 


time 


(a) Airfoil Unsteady - Freestream Airflow Steady 


i. Sinusoidal 
Vertical 
Translation 


ii. Horizontal 
Change in 
Motion 


iii. (a) Sinusoidal 
Change in 
Angle of 
Attack 


(b) Step 

Change in 
Angle of 
Attack 


y 

0 



time 



(b) Free Stream Airflow Unsteady - Airfoil Steady 


i. (a) Step 

Vertical 

Gust 


(b) Sinusoidal 
Vertical 
Gust 


ii . Horizontal 
Gust 


0 I 

\ > 

0 time 



COMPONENT OF LIFT DUE 
TO NON -CIRCULATORY 
(POTENTIAL) FLOW-INERTIA 


0 


1 

IT 


irpc 2 y 


COMPONENT OF LIFT DUE 
TO CIRCULATORY FLOW 
(VORTEX SHEDDING) 


DEFINITION OF TERMS 


-y 


c 



^ irc 2 p(Ud - ^AN) 


^ npc 2 UAa6(0) 


itpcU 2 F{ k) [a + - A)) 


irpcU 2 Aa^(T) 


irpcU 2 j *(^jp-) 





As indicated by the response functions in Table A-l, the change in lift 
characteristics is delayed in response to instantaneous change in airfoil angle 
of attack or flow velocity. The two physical mechanisms causing delay are 
found to be: 

1. inertia of the mass of air requiring adjustment and 

2. the finite time for the trailing edge vortices to be shed and cause 
adjustment with the flow field leading to a damping of the unsteadi- 
ness. 

These phenomena are analogous to a damped mass/spring system which can be 
described by a second order differential equation, such as: 


d 2 x 

m it T + 

dx 

c d£ + 

c 2 x - f(t) 

where: 

m = 

mass 


C 1 = 

damping factor 


C 2 = 

spring constant 


(A-10) 


An unsteady flow field about an airfoil will most likely be composed of a 
combination of several of the above motions. Nevertheless, each of the above 
motions can be approximated by a second order differential equation and the 
composite flow field modeled by a general second order equation. This response 
of the airfoil to these motions in unstalled flow will form the basis upon 
which to develop the phenomenological model of an isolated airfoil subjected 
to angle of attack excursions beyond the steady state stall limit, and inti- 
mately to represent the phenomenological model in usable, mathematical terms. 
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APPENDIX B 


Solution of the Differential Equation for the 
Effective Angle of Attack 


The effective angle is related to the instantaneous angle of attack 
through a second order linear differential equation. The differential equation 
is solved below for step, ramp, and sinusoidal changes in the instantaneous 
angle of attack. LaPlace transform techniques are utilized in the solutions. 

In terms of dimensional time, t, and dimensional time constants, t^ and tg, 
the governing equation is 

d^ (a - a ) + (~— + 7 —) 4r (a - a ) + r^T (a ec - a ) 

— =■ ' eff o '•t, dt eff o' t,t_ eff o 


dt 


12 


77 - (a. .. " a ) 

ins t o 


(B-l) 


It is convenient to non-dimens ionalize this equation by the airfoil velocity, 

U , and chord c. A non-dimensional time and time constants are defined below: 

T = tU/c 

T i ■ 'i u ^ 
t 2 ‘ C 2 U/c 

Employing these definitions Equation B-l becomes: 

“T <c, eff - “o> + (77 + 77) 37 - “o> + TV ( “eff ‘ °o> 

dx 12 12 


= — r 1 (a, ^ - a ) 

T 1 T 2 i nSt 0 


(B-2) 


This is identical to Equation 2 in the min text. Equation B-2 is a linear 
second order differential equation, is conveniently solved by LaPlace Trans- 
form techniques and is represented in the LaPlace domain as: 


L ( °eff * 

a o ) 


V 

where 

T (s) 

and 

s 


t s + 1^ t s + 1 ) T(s) 


(B-3) 


= LaPlace dummy variable 
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The effective angle of attack can be obtained for simple changes in the 
instantaneous angle of attack, Q!i ns t - & Q > by the method outlined below: 

“eff ' “o = L ’ 1 J L(a eff " “o^ = L_1 { ^“inst ~ a o ) x T(s) ^ (®-^) 


For example, assume that the instantaneous angle of attack, a inst> 
increases as a step function of time as shown in Figure B-l. 


1 ( “inst - “o ) ‘S' 


a ec - a = L 
eff o 


-1 


s (s t i + 1) (s t 2 + 1) 


a ,, - a 
eff o 


inst o 


n- - i + tv '' 7 ' 1 - V T,I » 


2 1 


(B-5) 


The effective angle lags the instantaneous angle as shown below. 



Figure 3-1 InposH of Effective Angle of Attack to a Step Change 


The actual shape of the a curve will depend on the system time constants, 
and r 2 > as indicated in (B-5). 

As another example the effective angle is found for a ramp input in a inst , 
with a ramp rate of increase, da = constant = Ct. 

dt 
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L <a inst ‘ “o ) ‘ c l /s ' 


a cc - a 
ef f o 


= L 


-1 


S (ST^ + 1) (ST 2 + 1) 


a .. a 
ef f - o 


= c 


1 T 2 - T 1 


-" T/T 2 - e ~ T/Tl + T - T 2 - \\ «-« 


An arbitrary variation in the instantaneous angle of attack can be repre- 
sented by a Fourier Series as discussed in the section "Extension to 
Arbitrary Variations in Angle of Attack." A Fourier series utilizes the sum of 
sine and/or cosine waves to approximate the variation. Solution of the differ- 
ential equation (Equation B- 2 ) for a sine wave will therefore allow solution 
for arbitrary instantaneous angles of attack. This solution follows: 

a inst” a O = ^ a inst sin ^ cot ) 


L (a. - a ) 

inst o 


a - a 
ef f o 


a - a 

ef f o 


A oi. 


inst 


" | 

(si^ + 1 ) (ST 2 + 1 ) (s 2 +(J 2 )) 



+ 




+ 4 k 



+ 



sin ( 2 kx - 4 > ) 




(B — 7 ) 


where k = 10/2 
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The exponential decaying terms are neglected for periodic motion, simplifying 
the amplitude of the sine wave results in 


eff o 

Aa . , 

in st 


sin (2kT — ) 

(1 + 4k 2 x 1 2 ) (1 + 4k 2 T 2 2 ) 


= f(k) sin (2kx - <j>) (B-8) 


The amplitude of the sine wave is noted as f(k) for convenience and shown in 
Figure B2. The phase shift is 


$ = tan 1 (2kx^) + tan 1 (2kT 2 > 


(B-9) 



Reduced Frequency, k -<jc/2U 


Figure B2. f(k) vs Reduced Frequency for Several Values of tj. 
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APPENDIX C 


Preceding page blank 


Increase in Airfoil Maximum Lift Coefficient for 
Unsteady Flow Test Data 


A limited literature survey of the effect of unsteady flow on the maximum 
lift coefficient of an airfoil was conducted to establish an estimate of the 
time constants associated with airfoil stall in unsteady flow. 

Available data were for a ramp rate of increase in angle of attack defined 
as an increase in angle of attack at a rate = d a /dt = constant. The lift 
coefficient achieved at stall with this type of unsteady motion was found to 
be greater than that achieved with steady flow. This increase, A is 

shown in Figure C-l as compiled from several sources. 

Solution to the governing differential equation, Equation 2, for the ramp 
rate of increase was given in Appendix B. The analysis is dependent on the non- 
dimensional time constants and T o . Based on a comparison between compressor 
component data and the analysis it was established that by setting the two time 
constants equal, r 1 = r 2 » the best match between data and analysis could be 
achieved. For this reason, analytical results are shown compared with the 
data in Figure C-l for values of = t« =t= 2.0 and 3.5. The comparison 
indicates time constants on the order of 3c/u are representative of the 
interaction between the unsteady flow and stall. 


PRECEDING PAGE BLANK NOT FILMED 
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STALL LUTE COEFFICIENT IM) 


symbol airfoil 


O NACA 2301X 

<1.3 

*'66,2-(l.8)(l5.2) 
° ">66,1-(1.8)(12) 

.6-13.4 

O Go 39o 

.36 

O Go 398,459 

.12-. 48 

V NACA 2R X 

•1-.3 

A - 

- 

O NASA 65,213 

1. 3-2.8 

A NASA 65,209 

2.8 


ANALYSIS 

r fV 3,5 
r V 2 ;V 


DIMENSIONLESS PITCH RATE, c d g-radians 

uat 

FIGURE C-l Increase. in Stall Lift Coefficient vs Dimensionless Pitch Rates for 
c M-rfYiil Tests 






APPENDIX D 


Relating Inlet Distortion to Rotor Blade Lift Coefficient 


The objective of this Appendix is to develop the relationship between the 
change in rotor blade lift coefficient required to perform the increased com- 
pressor work and the inlet distortion. Since a multistage compressor is built 
of individual stages haying similar characteristics, the effect of a change in 
maximum angle of attack of one stage will therefore be similar to the response 
of the next stage and derivatives or influence factors should be similar. The 
procedure will be to equate the change in work (due to distortion) done on the 
airflow to the change in work done by the rotor blade. From this equality will 
come the relationship between and pressure ratio which can be evaluated in 

terms of the operating point and the total pressure distortion. 

Change of Work Done on the Air . - The compressor work per unit time done 
on the airflow is: 

AE = w Ah 
a 


pA UAh 
x 


(D-l) 


where W a 
A h 
(AE) 
Ax 

U 


airflow = p 

ideal compressor total enthalpy increase 
increase in energy of the air per unit time, 
flow area normal to the axial direction 
axial velocity 


Consistent with the parallel compressor model which is based on a uniform 
compressor exit pressure, inlet total pressure distortion will require a change 
in work done on the air to maintain this uniform exit pressure. The change in 
work can be written: 


d (AE) 
AE' 


, dA 
d£. + x 


dU d(Ah ) 
U Ah 


(D-2) 


This change of work done on the air must be balanced by the change in 
work done by the rotor. 

Change of Work Done by the Rotor. - The rotor work is as follows: 


E RTR 2 PV R ArTR U RTR C L 


(D-3) 
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where Ep.pp = work per unit time done by the rotor 


Vr 


a rtr 


air velocity relative to the rotor 
rotor area = chord X unit span 
UpTp = compressor blade velocity 
Cp = rotor lift coefficient 

The useful work done on the air by the rotor will be somewhat less due 
to the finite efficiency. Thus: 


= nE T 


" 2 pV R A RTR U RTR C L n 


(D-4) 


The change in useful rotor work required by 

dp. 2dV R dA RTR dU RTR + + 

p V R ^StTR U RTR C L 

Equating the Change in Work . - Since the useful work done by the rotor 
must be equal to the work done on the air equations (D-2) and (D-5) can be 
equated: 


d( V 



distortion is therefore: 


dn 

n 


(D-5) 


dp . ^Sc , dU . dAh 
~ + A + U + Ah 

K X 


dn , 2dV R , dA RTR , dn RTR , dC L , dn 

p V R ^SlTR U RTR C L n 


Since the rotor area, the axial flow area, and speed remain constant, and 
assuming the change in efficiency due to distortion is strictly a secondary 
effect and therefore negligible, the equality reduces to the following: 


dU dAh 
U Ah 



(D-6) 


It remains to find a relationship between dCp, dU, and dVp. This can be 
established from the typical stage velocity diagram. 
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Velocity Diagram , - In the compressor mid stages, the stage stator exit 
angle will be to a first order independent of distortion around the circum- 
ference. Making this assumption the velocity into the following rotor will be 
as below: 


\ \ \ \ A A 


/ / 



Stage n-1 stator row 


Stage n rotor cascade 


T = a + $ • = constant 

tan 6 = tan (r - a) = U / U RTR 

U = U RTR tan (r - a) 

Taking the derivative of U 

dU = -U RTR (1 + tan 2 (T - a)) da 

and 

dU / 2“ a tan (T - a) \ dq (D-7) 

U ytan (T - a) 1 ) q 

Similarly: 

V R ' U RTR <COS <f - “»' 1 

dV R ’ " U RIR ( “ S (r " " 2 sln <r ' o) ) d « 

.tZS- . - ,»(!-.) ^ <»-« 

V a 
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Recall: 


C = 2ttci 

Li 


and 


dC T 


da 


Using these relationships equation (D-6) becomes: 


- a 


^ tan (f - a) 
C onib in ing t erms : 


- a tan (r - a) 


da + dAh _ 

a Ah 


= - 2 a tan (r - 


a) + 
a a 


da r 

a [ 


- a 


tan (r - a) 


+ a tan (r - a) 


-1 + = 0 
Ah 


Typical flow angles ( r — a ) are 35 degrees (Reference 13) and angle of 
attack of 12 degrees. 


tan 35° = .7 


a = 12° = .21 radians 


da 

a 


t 


.15-1 


1 


dAh 

Ah 


In essence 



dAh s da _ dC L (D-9) 

Ah a C 

L 

The additional required work comes from an increase in blade angle of attack 
(or lift coefficient). This additional work must be related to the inlet flow 
distortion. 


Compressor Work and Inlet Distortion . - The average work done by the 
compressor on the airflow remains constant. However, that section of the 
compressor operating in the region of low inlet total pressure must increase 
the pressure ratio to pump the flow to the uniform compressor exit pressure. 

It is assumed that this additional work is divided evenly among the compressor 
stages and the stage work, " A (Ah^)," can be related to the required increase 
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in overall pressure ratio as shown below: 


Ah = C (T - T ) 
1 P T i T i-1 


(D-10) 


and 


dAh j 

Ah~ 


where Ahu 
T Ti 
T Ti-i 


C p dTx i 


i— stage work 


total temperature at the exit of the 1 
compressor stage 


.th 


total temperature at the inlet to the 1 
stage 

(exit of ((i - l)th) stage) 


.th 


■i c P < v i vi ) 

The change in overall work is therefore given by Equation D-ll. 


E d(Ah. ) 
d (Ah) = i=l 1 
Ah n 

E Ah 
i=l 


n 

E 

i=l 


d(Ah ± ) 


Ah, 


Ah, 


n 

E 

1=1 


Ah, 


(D-ll) 


Since 


dAh. 

i , 

Ahi 

n d 


dC^ 
C L 
CL Ah. 


d (Ah) = i=l L 
Ah n 

E Ah 
i=l i 


Equation (D-ll) becomes: 


n Ah. dC /C 

1 Li Li 

n Ah. 


(D-12) 


and therefore 


d Ah 
Ah 


overall 



stage 


(D-13) 
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Thus it is shown that the overall increase in compressor work is related 
directly to the increase in airfoil lift coefficient of a single typical stage. 
By assu min g this additional required work is divided evenly among the com- 
pressor stages the required increased work can be equated to the required 
increased pressure ratio as developed below. The overall compressor work, 
assuming an efficiency of 100$, is given by: 


Ah = 


C 

P 



C 

P 



T t 
■ t 3 

•t t 


1 ) 


c 

p 




1 


Taking the derivative the following results 


d ( Ah) 

dR 

P 




Therefore 


d ( Ah) 
Ah 




(D-14) 


(D-15) 


where R 
P 


overall compressor pressure ratio = 



By combining equation (D-13) with this result the relationship between the 
increased pressure ratio and required blade lift coefficient is derived. 


dC L = dAh _ Y dR p (D-l6) 

C T Ah Y-l R 

L 1 - rT P 

P 

Realizing that the increased pressure ratio required of the compressor is 
the negative of the change in inlet total pressure, (distortion) 



(D-17) 


equation (D-l6) produces the relationship between the stage lift coefficient 
(dC^) and the distortion (dP T2 ) as given by Equation D-l8. 
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(D-18) 


da 

a 



Y-l 

1 - R < 1 ' Y)/Y 
P 



This result is combined with the effects of unsteady flow on the stalling 
lift coefficient to establish the effect of distortion on the compressor stall 
margin. 
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APPENDIX E 


Computation of the Loss in Conpressor Stall Margin 
Computer Program Description 


I. ABSTRACT 

A computer program has been written to calculate the loss in compressor 
stall margin due to inlet steady-state circumferential total pressure dis- 
tortion. The program is written in Fortran IV and is operational on the LTVAC 
CDC 3300 computer. The mechanics of the program are discussed below, and a 
program listing included. Input instructions and examples of the input and 
output are given. 

The program input is the circumferential total pressure profile which 
can be of arbitrary shape. A Fourier Series is then fit to this profile and 
the coefficients are modified by a transfer function to account for the delay 
in dynamic airfoil stall of the rotors. The loss in stall pressure ratio is 
then calculated as an output. 

II. FORMULATION 

A. Mathematical Description 

1. Distortion Amplitude 

The distortion amplitude is defined from the input total pressure 
field as the difference between the average and the minimum pressure divided 
by the average pressure (P^ - P t min)/Pf The distortion amplitude is used to 
normalize the calculated loss in stall margin. 

2. Fourier Series 

The technique of predicting stall pressure ratio loss for 
arbitrary distortion profiles utilizes Fourier series fit of the profile as 
indicated below. 

OO 00 

p^(0) = a Q + I cos (n0) + I b^ sin (n0) (E-l) 


PRECEDING PAGE BLANK NOT FILMED 
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where 


Pj. (0) = total pressure as a function of circumferential position, 
n 


= harmonic number 
. tt 


n 


P t 2tT ■ P T ( 0 ) de 

i 71 -TT ( f ,-?) 

= — / P T (0) cos (n0) d0 

— 77 TT 

b n = - / P T (0) sin (n0) d0 

TT i 

-TT 

The Fourier series can also be represented by sine terms. This 
will facilitate calculation of the loss in stall margin as will be shown later. 


Let 

a 

n 

b 

n 

c sin 
n 

c cos 
n 

(<f>n) 

(<t>n) 


2 

r 2 

c 

a 

t b 

n 

n 

n 

*n = 

tan ^ 

(a/b) 

n n 

Then 

a cos 

(n0) + 

b sin 


n ' ' n 


(n0) 



sin (n9 + 4>n) 


The series is thus written as 


P (0) » c 

— ^ = 1 - E — sin (n0 + t)>n) 

P T n=P P T (E-3) 

where the minus sign is the result of a change in the limits of 
integration from, — 7rto ir , to 0 to 2 7r . 
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3. Compressor Loss in Stall Pressure Ratio for Arbitrary Distortion 
Profiles 

The technique to establish the loss in compressor stall margin 
is fully described in the section of the main text entitled "Compressor 
Analysis." Briefly, loss in stall margin is defined with the aid of the fol- 
lowing figure which represents the process on a steady state compressor map. 



The "clean" stall margin is defined as 


SM 



R 

P 


0 


0 


(E-4) 
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The stall margin with distortion is defined as 

R - R 
PE PG 

SM Dist = R (E-5) 

P 0 

The loss in stall margin is thus the difference between these 

two equations : 


ASM = 



(e-6) 


This loss in stall margin is normalized by the distortion level 
as defined in Equation E7 


Dist = (P - P x )/P T 
min 


(E-7) 


resulting in a normalized loss in stall margin. Equation E-8. 

^SM 

(p - p )/p 
min 


(E-8) 


The technique to compute this normalized loss in stall margin 
due to an arbitrary inlet distortion profile is briefly as follows: 

(1) The input distortion profile is represented by a Fourier 

series. 

(2) The response of the compressor to this pattern is equal 
to the sum of the response to each of the components of the Fourier series. 

This is computed by the program. 

(3) The highest compressor stall pressure ratio, during a 
complete revolution of the rotor, is used to compute the loss in stall margin. 

B. Program Input 

1. The program is designed to compute the normalized loss in stall 
margin assuming a second order response with the non-dimensional time constants, 
r ^ = 7*2 = 3.5. This is the recommended mode of operation. In such a case 
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the only required input is the distortion pattern and the compressor reduced 
frequency, k. 


However, if it is desired to change the non-dimensional time 
constant and/or the "order of response" the program is designed to accept these 
as input items on Card 1 (Refer to Figure E-l). 

2. An arbitrary number of total pressure values are input into the 
program in order of increasing circumferential position up to a maximum of 144. 
The corresponding circumferential positions, which are arbitrary, are required. 
As an option it is possible to use the sine input subroutine to generate single 
or multiple harmonic sine waves of a specified amplitude in percent. Three 
sample cases are given which demonstrate the use of both types of input. 

Multiple cases are run by repeating all cards except the initial 
card which gives the number of cases to be run. 

3. The number of harmonics used in the Fourier series can be changed 
by changing the value of the parameter KN in the main program. KN = 72 is 
recommended for complex patterns where five degree resolution is adequate. If, 
for example, 2 l/2 degree resolution is required, set KN equal to 144. More 
harmonics will increase the accuracy of representing the input wave at the 
expense of more computing time; the effect on the calculated loss in stall 
margin is usually small. 

4. The input data arrangement is illustrated in Figure E-l. The data 
on the first card is used to signal the number of cases to be executed during 
the run. All control parameters are input for each case. The input data for 
the arbitrary distortion profile and the sine wave distortion profile are dif- 
ferent; and are illustrated in Figures E-l(b) and E-l(c), respectively. A 
detailed description of each input parameter follows. Parameters input in 

I format (see Figure E-l) are right adjusted. 


JJ 

NORD 

KREF 

T 1 

t 2 

LSFR 


Number of cases to be executed during run, input only once. 
Order of system response. 

Compressor reduced frequency. 

First time constant; if zero input, defaults to 3.5. 

Second time constant; if zero input, Tg defaults to 3.5. 

Key for long print out; LSPR / 2, long printout 

LSPR = 2, short printout 


KKKK - Key for sine or arbitrary distortion profile 
KKKK = 1, sine input follows 
KKKK / 1, arbitrary input follows 

K - Number of data points in arbitrary distortion profile defini- 

tion 1 < K < 149. 
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• • 


AND SO ON THROUGH CASE " JJ" . 


FORMAT 

PARAMETER 


FORMAT 

PARAMETER 


FORMAT 

PARAMETER 



/ FIO.5 1 F10.5 1 F10.5 

( T (1) T ( 2 ) 

' F10.5 1 

Tf 

FIO.5 

M 

! FIO.5 
T (K) 

F10.5 1 

FIO.5 

/ F10.5 1 F10.5 ! F10.5 

f PT (1) PT (2) " 

1 F10.5 1 

tt 

FIO.5 

TJ 

1 FIO.5 ' 
PT(K) 

F10. 5~ 

FIO.5 

l/l 6 ^ 
K 

/blank CARD 

/ 16 1 Fio.5 1 F10.5 ! 

NORD KREF t 1 

FIO.5 1 15 
t 2 LSPR 






Figure E-l(b). Arbitrary Distortion Profile Input Data Deck Arrangement. 


/ 16 1 FIO.5 
I mmm ampl 


/ 16 
/ KKKK 


|/I6 

NORD 


F10.5 

KREF 


T 


FIO.5 

T1 


Fio.5 r I5 

t 2 LSPR 


Figure E-l(c). Sine Wave Distortion Profile Input Data Deck Arrangement. 
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PT(I) - Total pressure distortion profile; arbitrary units, either 
ratio of local to freestream or avg, or absolute dimensions 
(PSIA, PSFA, N/M 2 , etc.). 

T(l) - Circumferential location (degrees) to describe distortion 
profile, must be input in ascending order. 


,MMM - Sine wave input distortion, number of cycles per revolution. 


AMPL - Sine wave distortion maximum amplitude. 


5. Most of the output data are labeled. The printed data include a 
summary of the input data, a short or long output data format, and a printer 
plot of the instantaneous and effective pressure ratios. 


SYS ORDR 
KC 

TAU]_ 

TAU 2 


For the short printout, the output data include: 

- System order, input parameter NORD 

- Compressor reduced frequency, input parameter KREF 

- System time constant one, T]_ 

- System time constant two, t 2 


SINE WAVE INPUT 


a. AMPLITUDE - Sine wave distortion profile maximum amplitude, 

input parameter AMPL 

b. MULTIPLE/REV NO. - Number of sine wave cycles per revolution, 

input parameter MMM 


or 

INPUT PT THETA - list of arbitrary distortion profile total pressure and 
angle pairs; these are PT(l) and T(l) input parameters, 
from I = 1 to K. 

PT AVG - Average total pressure of the input distortion profile. 

LOSS IN STALL MARGIN - This is the loss in stall margin. 

NORMALIZED LOSS - Loss in stall margin divided by distortion, where 

distortion is (Pr - Prmin)/ P T* 

The printer plot is composed of: 

1. the input instantaneous total pressure ratios 

2. and the effective total pressure ratios versus circumferential loca- 
tion, in degrees. 
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The long printout also includes 


PT/PTAVG 


THETA 

A(I)/AVE 

B ( I ) / AVE 


C ( I ) / A VE 


PHI (I) 

PHASE ANGLE 


- linearly interpolated total pressures divided by the 
average; computed from the input profile. 

Circumferential location, radians, 

- Fourier coefficients, a, for harmonics 1 to 18, divided 
by average total pressure. 

- Fourier coefficients, b, for harmonics 1 to 18, 
divided by average total pressure. 

- Fourier coefficients, c, for harmonics 1 to 18, 
divided by average total pressure. 

- See Equation (E-3) } radians 

- Transfer function phase shift for harmonics 1 to 18, 
radians 


AMPLITUDE RATIO - Transfer function amplitude ratio for harmonics 1 to 18 
PTp INST - Instantaneous total pressure computed from Fourier fit 

THETA - Circumferential angle, radians 


THETA - Circumferential angle, degrees 

PT 2 INST - Instantaneous total pressure ratio computed from Fourier 

fit and normalized by average total pressure. 

PT 2 EFF - The computed effective total pressure ratio, normalized 

by the average. 

SPR LOSS - The loss in stall pressure ratio. 

INPUT RPl/RPO - Defined by Equation (l8) in main text. 

RP(EFF)/RPO - See Equation (17) and Figure 8, in main text, for 

definition. 


C. Subroutine Descriptions 
1. Subroutine LINR 

Subroutine LINR accepts the input total pressure and angle data. 
The angles (and corresponding pressures) are input in increasing order and may 
be unevenly spaced. The subroutine linearly interpolates the data to genera e 
total pressures every five degrees (KN = 72). 
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2. Subroutine A1F 


Subroutine ALF accounts for the dynamics of the rotor stall pro- 
cess. The subroutine accepts (l) the number of harmonics to be used, (2) the 
order of the response, (3) the compressor reduced frequency, and (4) two time 
constants. The output is the phase angle and amplitude for each of the har- 
monics. 


3. Subroutine Sine 

The sine subroutine may be used to generate sinusoidal input data. 
When the sine subroutine is used, no pressure or angle data are input. The 
subroutine accepts the amplitude, and the number of sine waves per 3&0 degrees. 

4. Subroutine Plot 

The PLOT subroutine generates a printer plot of the interpolated 
input total pressure ratio, and the calculated effective total pressure ratio. 
The plot is suppressed when the output circumferential spacing is other than 
five degrees (i.e., KN * 72). 

D. Program Limitation 

The dimension statements allow for a maximum value of KN = 144 which 
corresponds to 36 harmonics, and l44 circumferential spaced total pressures. 

The plot limits are fixed; the upper limit is 1.25 and the lower limit is 0.75. 
If the order of response is other than two, the two input time constants have 
to be equal. 
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IV PROGRAM LISTING 


A listing of the main program and the four subroutines are given below: 
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PROGRAM 04815 TRACE 


COC 6600 FT N V4.0-P332 OPT=l 03/14/73 14.40.48. 




T 


"TD C 

A C" 

IT 

C 

15 C 

C 

c 

r 

c 

20 C 

c 

c 

c 


25 


30 


35 

C 


ur 


45 


_C. 


50 


PROGRAM 04815(IMPUT,0JTPUT,T APE5=INPUT,T APE6=OUTPUT) 

REAL K,N,KREF 

DIMENSION A<50), B(50) , CC50), PHI<50>, PT(15Q), T<150), TT(150) t 
^PT“Aa50l , PTC<150> , PPC15 0) , TTT<15T>, THA (150 ) ,RPI < 15 J ) *RPEF (153) 
DIMENSION PHE(50>, CO r K(50), C TCC<150>, DSM(150> 

READ <5,17t)> jj 
WRITE (6, 3201 “ 

DO 160 11=1 *JJ 

CiCCUL ATI ON OF NORMALIZED STALL PR RAT LOSS FOP A RHP I T R A RY TOTAL 

PRESSURE INPUTS 

CALCULATION OF EOURTER SERlE S FOR ARBITRARY PATTERNS 

INRUTH ANGLE- ONLY ONCE 

KN = NO O E INTERPOLATED PROBES 

PHE- PHASE ANGLE FROM TRANSFER FUNCTION 

°TCC= INSTANTANEOUS PT ON COMPRESSOR MAP , ANGLE = EFFECTIVE ANGLE 
INPUT IS INTERPOLATED TO 72 FARTS FOR 5 DEG RESOLUTION.. 
DETERMINATION OF MAX NO OF HA c MONICS AND MIN STEP SIZE 
36'0/< MAX NO DF HARMONICS) = TH-TA CYCLIC 
THETA CYCLIC/4. = MIN DELTA THETA. 

EX. 18 HARMONICS FOR MIN DEL T n THETA OF 5. OEG I.E. KN=72. 

JJ=NO OF CASES* TH= ANGLE, PT fl) =TOTAL PRESSURE ,K=NO OF RAKES 
NORO= ORDER OF SYSTEM RESFONS E *KREF=COM PRESSOR REOUCED FREQUENCY * T1 
AND T2 ARE TIME CONSTANTS* L^PR*? FOR SHORT PRINT OUT 
KN=72 
KNH=KN-f 1 
MM=KN/2 
KNN-KN 
QQ = KN 
EE = 0. 9 ~ 

RE AO ( 5 f 17 2) NORO , KREF, T1 , T? , L SPR 
IFdl.EQ.O.) Tl=3 • 5 
IFCT2.EO.0.) T2=3.5 
WRITE (6, 174) 

WRITE (6, 172) NORD*KREF *T1*T2*LSFR 
RE AO( 5 *17 01 NKKK 

KKKK=1 FOR SIN NTHETA , MMM= NO • /REV 
TF(KKKK-l) 4*2,4 
2 READ! 5,195)HMM,AMPL 

WRITE (6*196) KKKK,MMM, AMPL 
CALL SINE < AMFL *MMM,KN, FF, TTTl 
4 CONTINUE 

IF (KKKK-1) 8,35,8 

ft continue 

READ <5,180) K 

L = K 

KK = K 

INPUT FTRST ANGLE AT ZERO DEGREES 

WRTTE (6« 310) 

WRITE (6*19?) 

REAP <5*198) (PT ( I) * 1=1 » L) 

WRITE (6« 310) 

REPO (5,1 99 ) (TCI), 1=1 * L ) 

WRITE <6,190 (PT(I) ,T(T) ,I=1,L) 

nPTl=50. 

DO 20 14- 1 , L 


PAGE 1 
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PRO CRM D4B15 TRACE 


COC 560 0 FTN V4,0-P332 OPT = l Q3/14/?3 1 4 ,40.48* PAGE 2 


IF (PT<I4)-nPTl» 10,21 ,2- 

10 0PT1-PT (14) 

h ?q CONTINUE 

o qp 30 16=1, KK 

60 Tf 161 =T (16) /57. 2957795 

70 CONTINUE 

CALL LINR (PT,T,KK,KN, PP, TTT) 

75 CONTINUE 

DO 40 19=1, HH _ 

6 5 J1=?*I9-1 

J2=J1+1 

JT=J2+1 

EE = EE»PP(Jl)*4»fl*PP<J?) + PP(J3) 

4 Q- CONTINUE 

_7_r nTH-6.2B31853/QQ 

AVE-£E*nTH/6,2S31853/3 , 

WRITE (6,310) 

WRITE (6, 2 08) 

WRITE <6, 21") 

75 WRITE (6, 220) AVE 

WRITE (6*310) 

IF(KKKK.N E «1) GO TO 49 

DPT1-20Q OO P. 

DO 46 KI -1 , KN* 

8 0 IF (PP ( KI) -PPT1) 44,48,48 

__ 44 DPT1-PPCKI) ' __ 

48 CONTINUE 

49 CONTINUE 

P FT1=DFT1/AVE 

85 QPT1=1.Q-DPT1 

DO 50 LH=1,KNM 

PT(LH)=PP(LH)/AVE 

T(IH) =TTT(IM) 

: 50 CONTINUE 

90 I F ( L S PR* 2151*52.51 

51 CONTINUE 

WRITEtfc.igM 

W»IT€ (6,2001 <LP,PT(LP) ,T(LP) «LP=ltKNHI 

52 COMTINUE 

95 C IT TAKES 4. DATA POINTS FOR CORRECT INTEGRATION 

HMIMH/2 

DO 70 M=1,HHL 

N=H 

A A - 0 . Q 

_10Jj 89=r0.0 

DO 60 12= It 

J4=?»T2-1 

J5= J4+1 

J6=J5+1 

1 O F ANG1=N»T( J4) 

ANG?sN»T( J5) 

ANG3=N»T ( J6) 

AA = AA»PT< Jj »)»CQ5(flN .Gl> *4, »PT ( J5) »C0S. (ANS 2) »P TtJ6) T0 SfA NG 3) 

&fl=BB±P T 1 SIM LQMGl) t_4.. . »P T ,( J5).!tSl!LK A NS 2) »PT( J6) » S IN ( I> M11 

110 60 CONTINUE 










PPOGRAM 

D4815 

TRACE 

CDC 6600 FTN V4.0-P332 OPT=l 03/14/73 14.40.48. 

PAGE 

3 


C CALCULATION OF FOURIER COEFFICIENTS 

A« N):rAA*DTH/ 3. 14159265 36/3.0 
MN) = 0B*DTH/ 3. 14159265 36/3.0 
CTN^SaPT (A (H)*A(M)+B< M) *0<M) ) 

IIP D=A(H)/B(H) 

PHI(H)*ATAN(0) 

70 CONTINUE 
IFCLSPR-2) 71,72,71 

71 CONTINUE 

1 ?T WRITE (6,310) 

WRITE (6,212) 

WRITE 16,270) 

WRITE (6,240) (H, A(H ) , 0 ( M) , C (Mi'T'FHK M) ,N = 1,MML) 

72 CONTINUE. 

125 C CALCULATION OF PT(IV 

BCHAX=0. 0 
BC*AX=AHAX1(B> 

IFCBCHAX.EO.fl.O) BCMAX=-1.0 
SIGN=BCMAX/A1S(BCMAX) 

130 00 150 LK = 1 , 1 

IF(Tl-O.O) 76,74,76 
74 CONTINUE 

C T1=T2 IF N NOT =2 

IF1T2-0. 0)76*75, 76 
ri5 75 CONTINUE 

Tl=3. 5 
T2=3.5 
76 CONTINUE 

CALL ALF ( HHL , NOPD , KREF , Ti , T2 , P HE , C OF K) 

140 IF (LSPR-2) 77*78,77 

77 CONTINUE 

WRITE (6,250) 

WRITE (6,260) (LL,PHE< LL) ,COFK (LL) ,LL=1, NHL) 

78 CONTINUE 

145 C 18 HARMONICS ARE CALCULATED FOR KN-72 

THETA‘0.0 

00 80 KJJ=i,KNM 
THA(KJJ)=THETA 
THETA=THETA40TH 
150 80 CONTINUE 

C OUTPUT DATA ARE CALCULATED AT 5 OEG INCREMENTS 

DO 100 14 = 1 , KNM 

FTEF=9.0 

PTTxQ.O 

155 DO 90 15=1, NHL 

N=T5 

ANG=N*THACI4)+PHICI5) 

ANGE=ANG-PHE(I5) 

BA=C(I5)*SINTANG> 

160 PTT=PTT+B A 

C PTT IS THE LOCAL DISTORTION =CN»SIN(NTH + PHI) 

PTEFsPTEF+COFK( 15) »SIN(ANGE) »C(I5) 

QQ CONTINUE 

C ATE - PTT IF THETA 7E*0 =0.0 

M 165 PTCfl 4)=1 *0+ SIGN*PTT 

£ 





PROGRAM 

04815 TRACE CDC 5600 FTN V4.Q-P332 OPT=l 03/14/73 14. 40.48. 

PAGE 
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PTA(I4)=PTC(I4)*AVE 

TT (I4)=THA(I4)*57.2957795 



C 

PTCC(I4) IS THE EFFECTIVE PRESSURE 



PTCC( 14) =1* 0* SIGN*PTEF 

H 

17C 


R**I (I 4) =2 • 0 — PT (14) 



H 

ro_ 



REEF ( 14) =1.0 - SIGN*FT EF 



0SM(I4)=PTEF 



iro 

CONTINUE 



DSMAX=0.0 


175 


OO 120 14=1, KNM 



IF (OSH( I4> -OSMAX) 123,110,111 



110 

0SMAX=0SM(I4) 





120 

CONTINUE 



DSMN=0SMAX/0«T1 


1 nr 


IF (DSMN-1.0) 140,140,130 





130 

HSMN=1.0 





14Q 

CONTINUE 




WRITE(fc,31") 


WRITE (6,2701 

165 WRITE (6,26 0) OSMA"x,nSMN~ 

IMLSPR-2) 1*1,142,141 

141 CONTINUE 

WRITE (fe, 310) 

WRITE (6«?qPl _ 

190 WRITE (6,300 > 04 iP T M T 4 ) _, T H a ( 14 ) , TT (1 4 ) , P TC (I4),PTCC ( I_4 > , D SH ( 1 4 ) , 

114 = 1 , KNM ) __ 

14? CONTINUE _ ‘ 

WRITE (6,310) 

1E0 CONTINUE 

195 CALL PLOT ( KNH, RFI , RPFF , L SFP) 

WRITE (6,320) 

IE 0 CONTINUE 

STOP 

C 1 ~ 

2M 17Q format (IX, 15) 

172 rQPMAT(lX.I5.3Ein.5.T5 » 

174 FQPHAT(1X,60HS V S PROP KC T AU1 TAU2 LONG/SHOPT PR 

1INT OUT) 

1*0 FORMAT (IX. FI**. 5) 

?05 1°0 FORMA T C8F1P .S) 

1°? FORHATdX «1 5HINPUT PT THETA ) 

1°5 F0RHAT (1X, I5,F1H,5) 

1°6 FQPMAT(1X ,2T5,F1Q.5 ) 

19 6 FORMAT (9X . 1 6HPT /PT AVG THETA ) ~ 

210 2C_Q RQPHAT 1 1 X, 15. 2F1~ . 5 ) 

2 P 8 FORMAT (IX ♦ 11 HOUTPUT DATA) 

210 FORMAT (IX > 6HPT AVG) __ 

22_Q FQRHAT (1X.F1Q.5) 

230 FQRHAT ( 5 2HI A(I)/AVE e<I)/AVE C( I) /AVE PHUT) )_ _ 

215 232, FORMAT (1X.2QHF0URIER COEFFI C IENTS) 

24Q-.F0R HAT <1X,IS,4F10.5 ) 

2 r 0 FORMAT (4X.42HHARM0N IC NO P H *SE ANGLE AHPLITUOE RATIO) 

2EQ FORMAT(9X,I5.F10.5.5y. F I 0 .5) 

270 F0PWATC1X.4EHL0SS IN SIA_L L_MARGIN NORHALI 7E0 LOSS) 

220 260 FORMAT (1X.F1Q.5.20X. F1 0 . 5) " 





SUBROUTINE AIF 


TRACE 


COC 560 0 FT N V4.0-P332 Q°T=1 n i/ium 14,41.31 


PAGE 


1 




10 


15 


20 


75 


30 


35 


40 


SUBROUTINE ALE (HtNORJ ♦ KREF » TltTg.PHT.C 0 FK) 

C HsHARHONTC NO,COFK=F(K) 

REAL HAGO,HAS»K,KREF 

OIHENSION PHT ( 50 ) » CQFM 50 )* HAG (5 0 ) «_ 

C A AND B ARE RECIFOCALS OF NONOIMENSIONAl TINE CONSTANTS 

C A AND B DETERMINE THE STALL FILTER CHARACTERISTICS 

ORD=NORO 

Aal./Tl 

B«l./T 2 

WRITE ( 6« 1601 

WRITE < 6 t 1401 

WRITE < 6,1501 T 1 ,T 2 ,KREF 

IFIT 1 -T 2 ) 75 , 15,75 

15 CONTINUE • 

IF CNORO -21 17 * 25»17 

17 CONTINUE • 

OO 20 J= 1 1 M 

H*J 

k=kref»h 

ARGg 2 .»K»Tl 

PHI(J)=ORQ»ATAN< ARG) ; 

WAGQsfl»» 4 .»rH<»Tl»Tl) 

EXPsQRO /2 • 

HAG(J)= 1 . / f HAGD**EXP) 

COFKI J)=HAG<J) 

70 CONTINUE 

GO TO 135 : 

_ ?5 CO NTINUE 

DO 1 ? 0 J=l t H 

Hsj 

k=*ref»h 

Ay= 2 .»K/A 

A 4 = ?.»K/B 

PH HJ)=ATANt A 3 ) + AT AN (A 4 ) , 

HACP-SQ R T L.L l rO~ 4 **K*K/T A»B) I (A*B) ) ♦ 4 . 

1 <A»B»A»P) ) 

HAG(J)= 1 .Q/H AGO 

COFK( J) =H AG ( J) 

IT 0 CONTINUE , 

135 CONTINUE 

WRITE (fe . 16 01 

RETURN 

_C 

14 Q F 0 RHAT(*X. 2 QHTAU 1 TAU? KC) 

150 FORMAT ( 1 X. 3 F 10 . 5 ) 

lfO FORMAT (IX , n 

JlMJD 
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SURFOUTINE SINE 


TRACE 


con 6610 FT N V4 » 0 -P 3 32 QPT=1 03/14/73 14.41.07 


PAGE 
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SUPPQUTINE SINE(AHPL,>1MM,KN, PP.TTT) 

PI PENSION PPI15 01 »TTT(15Q) 

ft = KN 

KNM=KN+1 

DTH- 6.2831853/A 

PMWspMW 

WRITE (6<fe61 

WRITE (6t65> 

WRITE (6« 6 2) 

62 FORMAT (IX » 26HAMPLITU0E MULTI PL Ef REV NO.) 

WRITE(6«fe3) AHPL t PHH 

63 FQPHAT<lX,gF10»5) 

TH = 0. Q 

00 60 1 = 1 , KNH 

TTT<I1=TH 

ANGL=PHM»TH 

PPCI>= 1.0 - A WPL*SIN( ANGL) 

TH=TH+PTH 

60 CONTINUE 

65 FORMAT (IX i!5HSINE WA VE INPUT) 

66 FORMA T (IX «/) 

RETURN 

ENO 



SUBROUTINE PI OT TRACE CPC 6600 FT N V4.0-P332 OPT=l 03/14/73 14.41.10* PAGE 1 


SUBROUTINE PLOT ( KNM ,P T , RPE , LSPR) 

DIMENSION FTT150) ,RFE( IE 0 ) , N N (71 , 144) t LN Cl 20 ) 

QATA(NB=1H ) , <NI=1HI) , (NQ~1H E) 

IF(KNM-731 183,9,183 

_5 9 CONTINUE 

IF (ISPR-71 11,11,11 

11 CONTINUE 

WRITE (6,10) 

' 10 FORMAT (IX , 27HINPUT RPI/RPQ RP< EFF) /RPQ ) 

V> WRITE (6, 12) ( FT ( 1 1 ,RFE( I) ,1=1 , *NH> 

1? FQPHAmx,F10«5,10X,F10.51 

WRITE (6, 651 

13 CONTINUE 

PTH»1 • ?75 

_15 WRITE (6,781 

WRITE <6,161 

16 FORMAT (40X, 16HRPI/RPQ =INPUT=I) 

WRITE (6,171 

17 FOPHAT<4QX,70HRP<EFF17RPO=OUTPUT=E) 

_70 WRITE (6,781 

OO 70 J=1 , 21 

DO 18 I=1,KNH : 

NN(J,I) = 0 

18 CONTINUE 

?5 70 CONTINUE 

DO 40 J=l,?l 

PTM-FTH- » *)25 

PTl=PTH-0 #0 25 

DO 30 1=1, KNH . 

30 IF(FT(I) . GT ,FTH1 GO TO 75 

IF <PT (II , LE ,PTL) GO TO 75 

NN ( J, I) =1 

75 IF (RPE(I1 .GT .PTH) GO TO 29 

IF(RPE<I) .LE.PTL) GO TO 79 

36 NN ( J, I) =2 

GO TO 30 

79 CONTINUE 

IF (NN (J, 11-11 32.34,32 

72 NN(J,I)=0 

40 34 CONTINUE 

30 CONTINUE 

40 CONTINUE 

DO 76 J-l ,21 

DO 60 1=1. KNW 

46 IF (NN < J, H - 11 45 .50,55 

45 LN(I)=NB ; 

GO TO 60 

50 LN (II =NI 

GO TO 60 

_5T 55 L N ( I) =N O : 

6 0 CONTI NU E 

, n=.i 

Dl=0/1. 

02=0/11. 

55 D3=0/2l . 




V. SAMPLE PROBLEM 


A. Input Data Description 
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B. Output Description 

Pages 120 through 125 present the long output for the arbitrary distortion profile 
input. A two per revolution sine wave short output is presented on page 126. A 
final short output for a complex distortion profile input is presented on page 127. 
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SYS ORnR 

KC 

TAU1 

T AU2 

LONG/SHORT PRINT OUT 


2 

.10000 

3.50000 

3.50000 

0 



I NPUT PT THETA 


1_S 0 0 000 0 .90000 1. 00000 .90000 89.00000 1 .00000 90.00000 

1 . 2 0 000 91.00000 1.20000 359.00000 


OUTPUTOATA 

PT AVG 

1.12 407 



pt/ptavg 

THETA 


1 

.88962 

0 


2 

.80066 

P 08727 


3 

.80066 

P 17453 


4 

.80066 

P 261 80 


5 

.80066 

.34907 


6 

.80066 

.43633 


7 

.80066 

.52360 


8 

.80066 

.61087 


9 

.80066 

f 698 1 3 


10 

.80066 

.78540 


11 

.80066 

.87266 


12 

.80066 

.95993 


13 

.80066 

1 P 04720 


14 

.80066 

1.13446 


15 

.80066 

1.22173 


16 

.80066 

1 f 30900 


17 

.80066 

1 f 39626 


13 

.80066 

1 r 48353 


19 

.88962 

1 .57080 


20 

1.06755 

1 r 65806 


21 

1 .06755 

1 .74533 


?2 

1 . 0675=: 

1 .83260 



1.06755 

1 .91986 


24 

1 .06755 

2 1 007 1 3 


___ 25 __ 

i. 06755 

2 P 09440 


26 
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APPENDIX F 


Preceding page 



Solutions of the Navier Stokes Equations for Vortex Flow 


Definition of several vortex flow fields are developed below by solution 
of the Navier-Stokes Equations for both a steady and for a time dependent vor- 
tex. The choice of a vortex for use in describing turbulent flow will be made 
and the boundary conditions leading to this choice discussed. 


Navier Stokes Equations of Motion ; Using the cylindrical coordinate sys- 
tem shown below 



the equations of Motion of an incompressible viscous fluid in two dimensions 


are: 


3Vr . „ 9vr . v e 3Vr ve 2 . r _ I H . JlV + I IZt. + * _ 2 121 - V r) 

1t +Vr ^r r39 r ?r p 3r + V ll? + r 3r + P <30^ ^30 Vr ' 


32a. + v. 3tt + ». + 2flS* . . - i. |L + i§ijlL + I |2ft f \ + 2 

3t * 3r r 30 r r p 30 (3r z r 3r r z v 30 z 


(F-l) 


3Zr. _ v \l 

30 0'j 


(F-2) 


Continuity: 

If (rv r> + ff®- - 0 (F ' 3) 

For the vortex model, Fq and F r are zero and the pressure P and velocities 
v r and ve are independent of 0. Further assume that no source or sink occurs 
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within the area enclosing the vortex, in which case v r - 0. Equation (F-2) 
then reduces to the following: 


ilfi- 


3t 



(F-A) 


S olution of the Navier Stokes Equations . - There are several solutions of 
equation (F-A), each representing a different set of physical boundary condi- 
tions. Five of these solutions are discussed. 


For a steady flow v Q is independent of time and equation (F-A) reduces to 
the following: 


3 2 vfl . 1 _ 1 

3r* + r Sr r 7 



(F-5) 


The two solutions are: 

(a) For forced vortex (solid body rotation) : 

v e - V 

This is shown in Figure F-l. 

(b) For a Potential (free) vortex: 
v 0 - k 2 /r 

This is shown in Figure F-2. 


(F-6) 


(F— 7 ) 


A combination of equations F-6 and F-7 has been used to approximate a 
viscous flow model. Specifically, the following form results and is shown in 

Figure F-3. 


V0 


kr 

1+r 7 


(F-8) 


This satisfies the Navier-Stokes equations when r<<l or when r>>l; not, however, 
at r - 1 (see Figure F-A). Viscous effects cause non-recoverable losses. These 
losses dictate the vortex velocity and energy decrease with time. As a conse- 
quence a solution containing viscous dissipation must be time dependent and no 
steady state solution that fits the boundary conditions imposed by real flows 
can be found. 


For unsteady flow, three basic solutions have been found. The first is «* 8 
follows: 


\ 



(F-9) 
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Figure F-l. Velocity Ratio of a Steady Forced Vortex 
(Solid Body Rotation) - — Equation (F-6). 
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4 0 




Figure F-3. Approximation to a Vortex Having a Forced Rotation 

Near the Center and a Potential Motion at Large Radii 

Equation (F-8). 
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Figure F-4. Comparison of Che Approximation with the Two Steady Solutions 
of the Navier-Stokes Equations. 
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This represents e vortex that at zero time has a potential flow field. At 
t > 0, viscosity effects begin to cause the vortex to decay as shown in Figure 
F-5. The radius, r shown in Figure F-5 is that arbitrary radius having 
the velocity, vq, at time equal to zero. This solution is given, for example, 
in Reference l4 page 8l. 


A second solution exists for a vortex that is started impulsively at time 
t ■ 0, having a strength concentrated within a zero radius (a Delta function). 
At time t : - 0 the Influence of this line vortex spreads. The normalized shape 
of this vortex is shown in Figure F- 6 , and compared with the other vortex flow 
fields in Figure F-7. This particular solution is attributed to G. I. Taylor 
in Reference ( 15 ). - 


Its equation is given below: 


V0 


- B-p-e 


4vt 


(F-10) 


It will be assumed that this model most nearly represents the type of vortices 
in turbulent flow where the Influence of a formed vortex is at first limited but 
increases radially with time. Reasons for this selection are discussed below 
in Appendix G. 

A set of -additional solutions to the Navier-Stokes equations can be found 
by the technique of separation of variables. These solutions are as follows: 

V 0 ■ e ^ 1 Z (u) (F-ll) 

where: Z(u) ■ Jj(u), **( u )* (u) , H^(u) 

Jj(u) - Bessel Function of the first kind of order 1 

j_j(u) ■ Bessel Function of the first kind of order -1 

Y|(u) “ Bessel Function of the second kind of order 1 

H|(u) - Hankel Function of the first kind of order 1 

H*(u)' • Hankel Function of the second kind of order 1 

u 2 - (X 2 r 2 )/v 


These solutions however are oscillating with radius and as a result do not 
fit the boundary conditions of the problem at hand. 
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Figure F-5. Velocity Ratio for a Potential Vortex Allowed to Begin 
Viscous Decay at Time Zero. — Equation (F-9). 
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Figure F-6. Velocity Ratio of a Vortex Started Impulsively at Time Zero 
Equation (F-10). 
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Figure F-7. A Comparison of the Vortex Velocity Ratio for Three Flow 
Fields Satisfying the Three Different Boundary Conditions. 
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APPENDIX G 


Boundary Conditions for Vortex Model 

l 

For the purposes of developing a fluid dynamic model of turbulent flow 
certain boundary conditions must be imposed on the vortex description that meet 
the physical conditions associated with the turbulent eddies. These are assumed 
to be as below: 

(a) vq must be zero at a vortex radius of zero (viscous forces predominate). 

(b) v 0 must be zero at a vortex radius of infinity. 

(c) vq must be continuous between 0 <_r <_» and satisfy the equations of 

motion. 

(d) The zone of influence of an eddy or vortex must be small at first, as 
when it first forms, and grow with time; as opposed to an eddy that is 
fully established at all radii, then proceeds to decay, i.e. the trans- 
fer of momentum is outward with time. 

A summary of various vortex flows is given in Table G-l along with a graphical 
representation of the vortex and its respective ability to meet the boundary con- 
ditions. The Figure numbers in the Table refer to Appendix F. As shown, only the 
vortex formed by the impulsive start meets all the required boundary conditions. 

For this reason it was chosen as a typical vortex in turbulent flow and will be 
used in further development of the fluid dynamic model of turbulent flow. 

A word of caution, however, the vortex proposed assumes a laminar viscos- 
ity coefficient. Bear in mind that any given vortex may have smaller eddies 
forming in the core which in themselves may produce "turbulent" flow in which 
case decay would be governed by an eddy viscosity and occur an order of magni- 
tude faster. This turbulence could also cause modification of the velocity 
profile. Because of the assumption to use a laminar viscosity coefficient or 
specifically a viscosity coefficient independent of both radius and time, and 
because of the assumed impulse start of the vortex, the option must remain open 
to select other vortex velocity fields subject to experimental verification. 
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APPENDIX H 


Details of the Selected Vortex Flow Field 


A vortex flow field that satisfies the unsteady Navier-Stokes Equations 
and the proper boundary conditions was selected as representative of those 
eddies typical of inlet turbulent flow. This particular solution was first 
given by G. I. Taylor in 1918 and represents a vortex formed instantaneously 
in undisturb e<jl flow. The influence of this impulse begins to propagate out- 
ward at t = 0 . Table H-l is presented in summary giving the velocity, angular 
momentum, vorticity, circulation, and static pressure distributions of the 
vortex along with the rates of decay of each parameter. Derivation of the 
characteristics and graphical representations follow. 

Velocity . - The vortex flow field selected for detailed study is defined 
below: 

_ r 2 

v 0 " B ^2 e (H-l) 


Where: v^ 


velocity in angular direction 


B = constant 


r = radius 


t = time 

v = kinematic viscosity = fi/p 


The radius, a, at which the velocity is a maximum can be determined by setting 
the derivative^ dvg/ dr, equal to zero or 

r 2 

2Br 2 


#Vfl. „ 


3r 


B 

t 2 


iL 

4vt 


4vt 


* 0 


Thus: 



or 


(H-2) 

1. 

Note that a grows at a rate proportional to, t 2 . Normalizing the vortex radius 
by a, Equation H-l becomes: 


l4l 



■p- 

ro 


PROPERTIES OF A VORTEX 





'0 " FT? 


PT* B a e 


Normalizing Vn by the maximum value, the velocity ratio is; 


, -^) 2 -Dl (H- 

This is shown graphically (for a given time) in Figure H-l. 

Vortex Angular Momentum . - The equation of moment of momentum (angular 
momentum) is as follows: 

My - ($ x r)dm 

where the arrow denotes a vector quantity, "X" a cross (vector) product, and 
"dm" an elemental mass. In the cylindrical coordinate system 

$ - v 0 $ 

•+• ■* 
r - r r 

m - p x Volume 

dm ■ prd0dr (unit depth) 

where $ denotes the unit vector in the 0 direction 
r denotes the unit vector in the r direction 


Therefore: 


y.«v»2ir 

M » ■lj'o V « r ° rd6dr . li. 


To integrate let £ - r 2 

d£ - 2rdr 


ar 

r 3 . ( r 2) 3 / 2 _ ^ 3/2 
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Thus 




m . w * n ~ 

Mu - ifP72/ 5 e dt; 

M) 

M„, - TrpB(4v) 2 (H-6) 

Note: The angular momentum is Independent of time and hence Is conserved. 

Because is conserved it would be a logical term in which to express 
the constant "B". 

Thus: 


B 




irp(4v)‘ 


(H-7) 


Vorticity .- Vorticity is defined as: 
31 ■ I x v 




ft 





r 

4vt 


•! 

Br 2 ) 

2vt 3 | 


t 


(H-8) 


At a given time the maximum vorticity is found by setting the derivative 
of vorticity with respect to radius equal to zero. Thus: 


r 2 _ r 2 _ r 2 

an 2B 2r ” 4\>t Br _ ~ 4vt . Br 3 . 4vt _ r> 

17 FBt* ’rt 2 ' + 47 2 ? 1 * 0 


r @n - nMax " a 

■ea - lalMax ■ 0 (H ' 9> 

Substitute this radius in equation (H-8) the maximum vorticity is: 

2b 

l 0 lMax " t 7 (H-10) 

The max imum vorticity therefore decays proportional to 1/t 2 . 

Substituting the radius at which the vorticity is maximum (Equation H-9) 
into the expression for vorticity (equation H-8) and using this maximum vorti- 
city to normalize equation (H-10), the following becomes the normalized vorticity: 
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(H-ll) 


n 


r 2 

(*s(f) -l)e 


I^IMax 

and is shown plotted in Figure H-2. 

Circulation . - Vorticity is associated with a local elemental area. This 
vorticity can be different at each spatial location in the flow field. By multi- 
plying each elemental area by its associated vorticity and summing the results 
over the total area, the circulation Twill result. Thus: 


r - J~ • dX=^*^ • ds 
ds 


(H-12) 


V 5 


rde t 


• v 0 rd0 » B p 


T » 2irB ~2 e 


4vt 



(H-13) 


The circulation approaches zero at infinite radius indicating that the 
area associated with negative vorticity exactly balances the area weighted 
positive vorticity of Figure H-2. The maximum circulation can be found by 
setting the first derivative equal to zero. 


r 2 _ r 2 

ar 4xBr _ “ 4vt Br 3 ” 4vt =0 

ar t 2 vt 3 

t ■ constant 

1 _ ri. . 0 

t 2 4vt 3 

r - * / 4vt - /2~a (H " 14) 

Maximum circulation occurs at a radius a. 

Normalizing the circulation by the maximum value results in the following: 


l Max 


i,V - ft <7> 2 - 1) 


This is shown in Figure H-3. 


(H-15) 
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Figure H-2. Normalized (Vortex) Vorticity Distribution. 



Figure H-3. Normalized Vortex Circulation 


147 







Static Pressure . - By the assumptions made in developing the vortex, the 
momentum equation in the radial direction reduces to: 


v ft 2 _ 1 3P 

r p 9r 


(H-16) 


Integration yields: 

P o - P r " AP "jT dP 



^dr 



P o - 




(H-17) 


This static pressure difference can be normalized by the minimum static 
pressure, which occurs at r*0, to obtain the following dimensionless equation: 




(H-18) 


The normalized static pressure ratio and velocity ratio are shown compared in 
Figure H-4. 


Time of Origin of Vortices . - The vortex flow field defined herein began 
by an Impulse function at time t * 0. In reality this is only an approximation 
to such vortex motion. Assuming the vortex motion (however the details of its 
beginning) can be represented by equation (H-2) the virtual origin can be de- 
fined at such time as the vortex radius "a" is zero. Thus: 


a ■ /2vt 

ao‘ 

Conversely, given a vortex of radius a 0 the time of origin must have been to - ^ 
seconds earlier. 


The times involved with vortices having reasonable radii, i. e., 

a 0 is greater than .003 meters, are so great that viscous decay for vortices 
formed in the inlet is virtually non-existent. To illustrate, the velocity will 
decay to a value 1/2 of a given value by the following: 


V0 t-t? 


V0 


( 


t-ti 


tl ) 3 ' 2 

t2 J 


*2 “ 


v e 




V 0 t=t 2 


2/3 


- ( 2 ) 


2 / 3 . 


= 1.6t, 
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VORTEX VELOCITY DISTRIBUTION 



VORTEX STATIC PRESSURE DISTRIBUTION 



Figure H-4. Vortex Velocity and Static Pressure Distribution. 


AMPLITUDE RATIO 



Figure H-5. Race of Decay for VorCex Flow Properties. 
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Assume: a-^ - .3 meters 

t l = \/ 2v 

t 2 = 1.6 t ± (from Figure H-5) = 6250 seconds 

The residence time of this vortex in a 3 meters long inlet would be =.02 seconds. 
For this reason the assumption will be made that the vortex is essentially in- 
dependent of time and defined as below: 

- 1 ( 1)2 

V 6 B" (£) e 2 a 

where: B" » B' /t 3 / 2 

B' - MvTv / (4p(4v) 2 ) 
a ■ v7vF 
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APPENDIX I 


Preceding page blank] 

Total Pressure and Flow Angle of a Vortex Superimposed on a Local Flow 


As a first step in developing the statistical model of turbulent flow the 
flow field of an isolated vortex was defined in the cylindrical coordinate 
system in Appendix F. In part A of the following discussion the vortex flow 
field will be converted to cartesian coordinates with the origin at the vortex 
center. In part B, the vortex is superimposed on axial flow in a channel and 
the flow properties of this moving vortex as measured by a probe fixed with the 
wall are determined by a transformation of coordinate systems. 

A. Single Vortex in Cartesian Coordinate System . - The following is a 
description of the vortex flow field velocity components, flow angle, static 
pressure variation, and total pressure variation as seen in the cartesian 
coordinate system fixed at the vortex center. This coordinate system is 
sketched in Figure 1-1, The circumferential angle, 9, is measured counter- 
clockwise from the positive x-axis. The vortex size, a, is the radius at which 
the maximum tangential velocity, vgoax occurs. These and other basic vortex 
relations to be used herein come from Appendix H. 


Velocity Components - 


by 


The circumferential velocity is a function of the radius only and is given 


Br -r /4vt 

7 * 


(i-i) 


where B is a constant dependent on the vortex strength (circulation) 


V is the kinematic viscosity 

t is the time of origin, but is to be assumed constant because of the 
short period that the vortex is in the field of interest. 


The maximum velocity occurs at the radius r 


v 0 max 





a = Jlvt 


( 1 - 2 ) 


PRECEDING PAGE BLANK NOT FILMED 
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SYMBOLS: 

a Vortex size,. v 0 = v 0max at r = a 

6 Circumferential angle, degrees 

x/a Horizontal non-dimensional coordinate 
y/a Vertical non-dimensional coordinate 

v Vortex tangential velocity 

u® Vortex horizontal velocity component 

v Vortex vertical velocity component 

o Flow angle 

r Radius from vortex center to flow field point 


Figure 1-1. Single Vortex Model in Cartesian Coordinate System. 
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So that in terms of the maximum velocity 



From the geometry of Figure 1-1, 
are 


„ ,r. - 1/2 [(r/a) 2 - 1] (1-3) 

the velocity components in the x and y directions 


u » -v 0 sin 6 

v ■ v 0 cos 0 


(1-4) 


The relation between the radius and the cartesian coordinates can be expressed as 

- 6 ■ 6 - fflt (I - 5 > 

2 

Substituting equations (1-3) 6 (1-5) into (1-4) and using the relation (r/a) ■ 

(x/a) 2 + (y/a) 2 yields 


“ ■ - V 9b« (y/a>e 


- 1/2 [ + (y/a) 2 - 1 ] 


(I-6a) 


v " V 6max (x/a)e 


- 1/2 [(x/a) 2 + (y/a) 2 -1] 


(I-6b) 


Flow Angle: The local flow angle is always tangential. This is due to the 

fact that no radial velocity exists in the proposed vortex model. The flow angle 
as determined from the vertical and horizontal velocity components is: 



“* 9 1 


’ -1 

a ■ arctan (v/u) ■ arctan 

■ arctan 

sin 6 

-'W 9 


cos 0, 


(1-7) 


arctan 


-1 

tan 8 


Static Pressure Variation: The absolute value of static pressure is required 

in order to later evaluate the total pressure. From Appendix H, the static pres- 
sure is 


4 - r 2 /2vt 
6 


(1-8) 
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where B, t, v have been defined, see equation (1-1). 

p la density . 

p is the uniform "freestream" static pressure at infinity 
P° is the static pressure at radius r from the vortex center 
r 

The maximum static pressure depression occurs at r ■ o, so 


«, - v- - 

This equation can be written in the form 


(P r - (W.V 


(1-9) 


By rearranging and squaring the expression for the maximum velocity, equation 
(1-2) yields 


(Ba/t 2 ) 2 - v 2 max e 


( 1 - 10 ) 


Substituting equation (I- 10) into equation (1-9) , together with the fact that 


vt ■ a /2 f yields 


^ P r ~ P o^max 


. v 2 

2 0max 


( 1 - 11 ) 


The ratio of static pressure at any radius to the maximum static pressure 
depression is 


_2 -r 2 /2vt 

>B v e 


P r“ P o 
< P r - P o>max 


-(r/a)‘ 


( 1 - 12 ) 


-PB Z v 

t 3 


as 


given in Appendix H. Therefore, the absolute pressure at any radius is 

-(r/a) 2 [1 “ (r/a) 2 ] 


p-p ■ (P - P ) e 
r o r o max 


- fi. e 

e 6max 


(1-13) 


Which, in terms of the cartesian coordinate system becomes 

2 [1 - (x/a)2 - (y/a)2] 

P r " P o " 2 V 0max 6 


(1-14) 
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Vortex Total Pressure: The total pressure in the vortex flow field, as seen 

in the vortex coordinate system, is equal to he sum of the static pressure and the 
dynamic pressure. The total pressure is 

P I ■ p r + f v e (I - 15) 

where 

P is defined by equation (1-14) and 
v n is defined by equation (1-3) 

D 

So that the total pressure becomes 


2 [l-(r/a) 2 J ? » [l-(r/a) 2 ] 

P o-i V Lax e + t V Lax< r/a > e , 

p + £ rr / ^ + t / ^ ii 2 U ^ ' (y/a> 1 ^- 16) 

p I + <y/«) -U v. _ e 


Velocity and pressure variations are shown in Figures 1-2 and 1-3. The results 
are plotted versus the dimensionless coordinate, X/a, for y/a = 0. The velocity 
distribution is maximum at x = a, is zero at x = 0, and diminishes to zero at large 
radii. The pressure distributions, both static and total, are shown as the ratio 
to the static pressure at infinity. The maximum static pressure depression occurs 
at r = 0. The total pressure reaches its maximum value at x/a s 1.4, approaches the 
local static pressure at r = 0 (V = 0), and approaches the static pressure at 
infinity at large radii (V = 0). 


B. Single Vortex in a "Fixed" Coordinate System . - The vortex flow field 
will now be defined in terms of a "fixed" coordinate system. The vortex model 
(and its coordinate system) are assumed to move at a constant axial velocity 
relative to a fixed channel as sketched below: 



157 
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4 


x/a 


tex Velocity Distribution 




P T /P 


o 


In the vortex model coordinate system, the flow properties are the same as 
previously defined. In the "fixed" system, however, the uniform constant velo- 
city must be added to the vortex velocity components. This yields different 
resultant velocity vectors and total pressures. In both cases, the static 
pressure at a point must be identical. (Properties of a flow field are indepen- 
dent of the coordinate system) . 

The fixed coordinate system is sketched in Figure 1-4. The vortex and 
its coordinate system are assumed to move at constant velocity, U 0 , which is 
defined as being parallel to the X-axis. The vortex system axes can be set 
parallel to the "fixed" system axes because the vortex flow field is dependent 
on radius only. The vortex flow field point of interest is the point which 
coincides with the "fixed" system origin. Therefore, the flow properties at 
the origin of the fixed coordinate system due to a vortex located at X, Y, can 
be obtained by applying the following transformation to the equations in the 
moving coordinate system, 

X ■ -x and Y ■ -y (1-17) 

The vertical coordinate remains constant, but the horizontal coordinate is a 
direct function of time 

AX - U 0 At (1-18) 

Velocity Components. - The vortex flow field velocity components are deter- 
mined by substituting (1-17) into (1-8), so that 


(I) . j s[(x/«) 2 +(y/») 2 -i) 
0 max V 


(I-19(a)) 


v 


,, (X/a) e -1 * f (X/a) z +(Y/a) 2 -ll 


(I-19(b)) 


And then adding the constant flow velocity, U Q 
U - U 0 + u 




(I-20(a)) 


V 


(Y/ a \ . -%[(X/a) 2 +(Y/a) 2 -l] 

< X/S > e 


(I-20(b)) 


Flow Angle. - The flow angle at the "fixed" system origin is defined as the 
arc tangent of the ratio of vertical to horizontal velocity. It is measured 
counterclockwise from the positive X - axis. The expression is 
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SYMBOLS : 

U q Flow field constant velocity 

a Vortex size, v Q = v flTnpy at r * a 

9 Circumferential angle in vortex system 

x/a Horizontal coordinate in vortex system 

y/a Vertical coordinate in vortex system 

v Tangential velocity 

u Vortex horizontal velocity component 

v Vortex vertical velocity component 

a Flow angle 

r Radius from vortex center to fixed system center 
X/a Fixed system horizontal coordinate 

Y/a Fixed system vertical coordinate 


Figure I- 1 *. Vortex Model in Fixed Coordinate System. 
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oi = arctan (V/U) 


a 


arctan 


-v e (X/a) e 
°max ' 


-h 


L U ° + V0 max (Y/a) 6 

which after some rearranging 


[(X/a) 2 +(Y/a) 2 -l] 

- 1 s[(X/a) 2 +(Y/a) 2 -l] 


a = arctan 


(Y/a)+Ua 


V0 


max 


- (X/a) 

e -4 [(x/a) 2 +(Y/a) i! -l] 


Static Pressure. - The static pressure in the "fixed" system is identical 
to the static pressure in the vortex system at a given point. Therefore, in 
the "fixed" system, the static pressure is 


Pr 


Po 



V0. 


max 


e [l - (X/a) 2 -(Y/a) 2 ] 


( 1 - 22 ) 


Total Pressure in Fixed System, - The measured total pressure in the fixed 
system is the sum of the local static pressure and a corrected dynamic pressure. 
The correction is due to the fact that the local (at the origin) flow angle io 
not aligned to the "fixed" total pressure probe. The resulting relation is 

PT * Pr + hR(a)(p/2)W 2 (1-23) 

where Pr is static pressure, equation 1-22, 


W is the resultant velocity vector, 
and nR(a)is a recovery factor which is a function of the local flow angle 
The resultant velocity vector is 

w 2 - u 2 + v 2 


W 2 - (X/a) 2 +(Y/a) 2 ve 2 ax 


- [(X/ a) 2 +(Y/ a) 2 -l] 


+ U G j U 0 +2ve max (Y/a) e f(X/a) 2 +(Y/a) 2 -l] j 


Substituting (1-22) and (1-24) in (1-23) and simplifying yields 

P T - P 0 + P/2 v 2 x e U“(X/a) 2 -(Y/a) 2 ]( r((x) [(x/a) 2 +(Y/a) 2]_ 1 j 




, x /9 „ i „ , v/a , ~h [(Y/a) 2 +(X/a) 2 -l] ) 

+ n R (a) p/2 u 0) u 0 +2v e max ( y / a ) e ( 


(1-24) 


(1-25) 

l6l 



Velocity components, flow angle, and pressure distributions are presented 
in Figures 1-5 through 1-9. For this single vortex in a uniform constant 
velocity flow, the dynamic pressure recovery factor, ) , is assumed to vary 

as the cosine squared function. The dynamic pressure recovery factor T )p( a ) 
corrects for the probe characteristics which at angle of attack yield total 
pressure lower than the actual. This variation with angle of attack is illus- 
trated in Figure 1-10. For this example, the vortex is assumed to rotate 
counterclockwise. The ratio of constant velocity, U Q , to the maximum vortex 
tangential velocity, V ftrn ny , is 1.2. The horizontal and vertical velocity 
component^ of the resultant velocity vector (vortex flow field superimposed on 
the constant flow velocity) are given in Figures 1-5 and 1-6. The coordinates 
(both X/a and Y/a) represent the distance from the "fixed" coordinate system 
origin, the sensing total pressure probe location to the center of the vortex 
flow field. The angle between the positive X-axis and the resultant velocity 
vector is given in Figure 1-7. These flow angles follow directly from the 
velocity components and determine the dynamic pressure recovery factor. The 
static pressure distributions is presented as Figure 1-8. It is basically 
dependent only on the distance from the vortex center. The total pressure 
distributions is presented in Figure 1-9. Note that whereas the velocity com- 
ponent distribution are symmetric for negative and positive Y/a, the total 
pressure distribution is displaced downward. This is due to the static pressure 
distribution and to the angle of attack recovery factor being less than unity. 
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Figure 1-6. Vortex Vertical Velocity Component in Fixed System. 
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Figure 1-8. Vortex Static Pressure Distribution in Fixed System. 



istribution in 
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FLOW ANGLE - a, Deg. 


Figure 1-10. Tvpica] Dynamic Pressure Recovery Factor , 
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APPENDIX J 


The Autocorrelation Function of a Random Signal' 
Composed of Several Independent Random Variables 


The autocorrelation function resulting from a stochastic process is 
found in functional form by statistical methods below. These developments 
are then applied to the vortex flow field and the autocorrelation function 
of the turbulence established. 

Autocorrelation Function of a Stochastic Process . - The objective of 
the following treatise is to establish the autocorrelation function of a 
resultant signal, f(t), composed of wave forms g(t) that occur randomly with 
time at an average rate of N per unit time (Poisson waves). Each waveform 
is specified by several variables that are random. Schematically; 



where g(t) may be structured as below: 
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The waves f(a, v, y, n) are identical for identical values of a, v, y and n. 
The resultant signal, assuming the variables occur randomly, will be: 



As a first step toward the development, assume the wave occurs randomly 
with time at a mean rate of n/second and has the waveform specified by only 
a single random variable, "a". 



Where a^ is governed by its probability density function "P(a). This wave can 
also be considered as composed of various sets of identical waves, such as 


Set 1 



Set 2 



where the number of pulses n^ of size a., and nj of size a p , etc. is established 
by the probability density function P (i^) p (aj). 
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The autocorrelation function of pulses of size a,, governed by a Poisson 
distribution having n^ pulses per unit time is (see p 336-7 of Reference (36)) 


R (t) 


n , R ' (r) +1 

1 a l 


(v^ 


(J-2) 


where R a ^ ( t) is the autocorrelation function of a single wave. 

Similarly, the autocorrelation function of pulses of size a ^ will be: 

2 C 


R (T) = n R ( 

a 2 Z a 2 


T> + ( f a 2 (t> ) 


(J-3) 


Continuing the process, the overall autocorrelation function can be obtained 
by summing the above autocorrelations. 


R a (T) " 2> a '“i < T > + lE\' a 

i 1 i 


(t)] 


■N/ P(a) Ra (x) da + l J ? ( a) f a (t) | da 3 (J_4) 


Subtracting the D.C. component (which is independent of (r)) from the 
autocorrelation function the following results. 


R a (T) 



P (a) R (t) da 
a 


(J-5) 


Extension to multi-variant signals 


The typical wave is represented by the following general function 
f(a, v, y, n) where a, v, y and n are independent random variables governed 
by their respective probability density functions: 

P(a) 

p(v) .... : . .. . . .. 

P(y) 

p(n) 
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Assume first v, y, and n are constant and only "a" varies. The resultant 
autocorrelation will be the same as given by Equation J- 5 . 


R a (t) = (7 ^ (x) + constant} v J 

i y \ constant 

n ) 


( J-6) 


or schematically 





constant 


Now, for a given a^, v will take on a random value illustrated schematically 
below: 




constant 


By Equation J-4 the autocorrelation of this signal of width a, and 
having random widths v^ as governed by a probability density will be: 


'v (T) -■£« R v (T) I 

a=a 1 a=a 1 


(J-7) 
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There are number of occurrences of this signal. Similarly for a 
■ i 


R v ( T > = £ “ 


i 

a*a r 


V i V i 


CO 


+ constant 


a=a r 


a 2 : 

(J-8) 


with n ' occurrences of this signal. Continuing for a - ag, 
. . ®2 a n the resultant autocorrelation function will be: 


* .V (,) 


®. 2" *„ < T > i (,) 


a l V i " V i 


a 2 V i " V i 


a=a. 


a=a„ 


+ . . . . n n R v ( T ) 


a i" v i " V i 


+ constant 


a=a. 


* v± «> 


+ constant 


(J-9) 


a=a 4 


Continuing this same process allowing y and n to take on random values, the 
autocorrelation of the resultant signal will be: 


‘ vayn 


< T > - E" y E»„E".lv "v, (T) + 


constant 


(J-10) 


For a continuous range of variables this will result in an integral form: 

r (t) = n/* f f f R CO P (a) P(v) P(y) P (n) dadvdydn (J-ll) 

va y n ••n* , y* , v* , a 


+ constant 

The constant is simply the mean of the resultant signal squared. Thus 

constant = [f (t)] 
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If these positive pulses are added to a wave train of similar but nega- 
tive pulses occurring at the same rate 


avyn 


(t) = 2N fjff R (t) P(a) p (v) P (y) P (n) dadvdydn 

+ f(t) (J-12) 


But since these negative pulses are equal in number and shape to the 
positive pulses and differ only in sign, f ft) = 0. Then the following auto- 
correlation results for equally likely positive and negative pulses of total 
rate N 


R 


avyn 


CO 


(T) 


P (a) P (v) p (y) p (n) dadvdydn 


(J-13) 


Autocorrelation Function of the Total Pressures . - The model of turbulent 
inlet flow is hypothesized as being composed of a random distribution of 
vortices each having a specific size, strength, direction of rotation, and 
location. The total pressure fluctuation created by each vortex is given by 
Equation 39 of main text. 


££x. 

<lo 


2 n(l£sS) 


*2 f(U 0 1/ a) ^(Y/a) 2 -l] 

(Y/a) e 




KU 0 t/a) 2 +(Y/a) 2 -l] 


( 39 ) 


For a specific vortex having a given set of properties (size, a; strength, 
v ©max> s P* n direction, n; and location, Y), Equation J-l4 signifies a single 
time function. However, each vortex has a different set of properties. There- 
fore, the flow field is composed of a family of time functions. The auto- 
correlation function of the total pressure fluctuation composed of the random 
vortices flowing downstream with the flux of N per second is given by Equation 
J -13 'withP(a), p(v), P(y), and p(n) being the probability density functions 
of the respective independent random variables, the autocorrelation function 
of the general wave R( ) is found by means of the definition of autocorrelation 
function for discrete waves. 


In Equation J-13, P(a), p(v), P (y) and p (n) are the probability density 
functions of the respective random variables as given by Equations J-lU through 
J-17. These density functions are in general described by a Beta probability 
density function. 
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This density function and its transformation are flexible and can be made to 
fit the boundary conditions of the respective random variables. See Appendix 
K for a detailed discussion of probability density functions. The respective 
vortex probability density functions are as follows: 


a. Vortex Core size, a 


'’(a) - f <f> 8 W 


n 


h' 


for 0 < a < H 


(J-14) 


b. Vortex Strength (Maximum tangential velocity), 

m 


, v Q m v n n 

k 0 a 6 v 

v 0 ) u V U J V u ' 

max o o o 


for 0 < v* 

o 


max 


< U o (J-15) 


c. Vortex Lateral Location, Y 


P (Y) 


1 

H 


for -h < Y < H - h 


(J-16) 


d. Vortex Spin Direction, n 

P (n) = +1,-1 (J-17) 

The autocorrelation function of the vortex flow field as measured at the 
total pressure probe is found by incorporating the total pressure wave 
(Equation 39 ), the definition of the autocorrelation for a descrete wave 
(Equation 4l), and the density functions into Equation J-13. This is given 
by Equation J-l8. 
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(J-18) 


This can be reduced as far as Equation J-19 in closed form. Integration 
with respect to the random variable f a, must be done by numerical techniques. 

A computer program was written to evaluate this integral. 
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where 

ERF(z ) 



Numerical evaluation of this integral yields the normalized auto- 
correlation function as shown in Figure 27 of the main text. 
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APPENDIX K 


Probability Density Functions 


A physical interpretation and definition of the probability density function 
”f(x)" will be given below using the common Gaussian density function as an 
example. A summary of the more common density functions and their characteris- 
tics will then be given. 

Transformation of the Beta distribution demonstrating the added flexibility 
that can be obtained by use of this technique will follow. Such flexibility 
is required to meet the wide range of density functions needed to fit the 
physical turbulence characteristics. 

Probability Density Function .- The probability density function "f(x)” for 
random data describes the probability that the data will assume a value within 
some defined range for a single event or, as below, at any instant of time. 
Consider the following time history of the signal x(t) below: 



E 

i » 1 


At. 


Prob [x < x(t) S x ♦ Ax] » lim 

T— oo 



The probability that x(t) assumes a value within the range between x and 
(x +Ax) may be obtained by taking the ratio of T x /T when T x is the total 
amount of time that x(t) falls Inside the range (x, x + Ax) during an observation 
time T. This ratio will approach an exact probability description as T approach- 
es infinity, often random data of this nature assumes a probability density 
approaching the bell shaped or Gaussian form as shown to the right of the above 
sketch. The mean value, not random data will be the average over the entire 
range. Thus: 

u.l. 

n * J' x f(x) dx (K-2) 

1 . 1 . 
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The mean square value, a , is the squared value of x over the entire range. 
Thus : 


u 

-/ 


u.l. 

r 

1 . 1 . 


x 2 f(x) 


dx 


(K-3) 


The mean and mean square values are called the first and second moments 
respectively. In general the m^h order moment is as follows: 

u.l, 

£ m = f X“f(x) dx (K-U) 

^ 1 . 1 . 


A random variable is completely specified by definition of all its moments. 

Examples of Probability Density Functions .- A summary of density functions 
is given in Table I along with a graphical representation. A typical source for 
these density functions is Reference 16. The characteristics of the functions 
are also given in the Table and will be significant when subsequently fitting the 
density function to the properties of the random vortex model. 

It should be noted that all of these density functions are greatly limited 
in flexibility. Specifically the moments, and hence the shape of the respective 
density curves, are fixed (with the exception of the Beta Density). 

As an example of experimental data consider the probability density func- 
tion of the total pressure fluctuations measured in an inlet of Reference 12 and 
as shown in Figure K-l(a). The data obtained from test agrees closely with the 
laus si an Density. However , at other conditions and measurement locations a 
skewed density has been measured in the same test. This is shown in Figure 
K-l(b). As evident, the Gaussian Density is not representative of this data. 

The moments of the Beta Density can be changed by choice of the constants 
b and c (see Table K-l). The limit values of this function remain, however, at 
zero and unity. The flexibility of this density can be greatly increased by 
transformation as developed in the following section. 

Transformation of the Beta Density Function .- The Beta Density Function is 
defined as: 


f(x) = Ax t (l-x) C for o g x g 1 (K-5) 

where A = r (b+c+2) 

T(b+1) T(c+1) 

and T = Gamma Function 
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Figure K-l. Comparison of Probability Density Functions Obtained 
from Inlet Turbulence Measurements with the Gaussian 
Density Function. (Test Data from Reference 12, 
page 321 and 326 respectively) 
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The Gamma function for integers is defined as 

r(n) - (n-1)’ 

Therefore for integer values of b and c the mean .value of f(x) is 

b+1 

f(x) " b+c+2 

For b = c = 2 the density function is as below: 



This can be skewed by choosing b / c and for example with b - 2 and c 3 the 
following density function results: 



To increase the flexibility of this density function it will be transformed in 
two steps. First f(x) will be transformed linearly by y = dx. By use of the 
following transformation (see for example Reference lb, r ) • 


f(y) = 


M 


dx 


(K-6) 


x d 


d 

The density function of y now becomes : 


f(y) = 




(K-7) 


l8l 



This function is shown below: 
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The first moment (mean) of this distribution is by -reason of the translation, 
z = y - y, equal to zero. The second moment, which is the mean square value is 
found by integration to be: 

o 2 - z 2 = d 2 h (b * 2 i 4^-r - riiiil] 2 ! (K-9) 

|(b+c+3J (b+c+2) L b+c+2 j 

Normalizing the transformed density function by the root of the mean square 
value, ct , results in equation (K-10) 

cf(z) = 4<ff + f) b <1 - <f f f + d» C ( K - 10 > 

To illustrate the flexibility of the transfermed Beta Density, equation (K-10) 
is shown graphically in Figure K-2 for various values of b and c, where b f c . 
Note the varying amounts of "skewness" that can be obtained. In Figure K-3, b 
is assumed equal to c which is equal to "n". Thus 

b = c = n 

This gives a symmetrical Density Function. For n = o, the uniform density is 
obtained; for n = the Gaussian Density results, demonstrating the wide 
range of density functions that can be formed from this transformed Beta Density 
Function. 

The transformed Beta density still has only the two constants b and c as 
variable to change the general shape of the density. As a result this density 
still falls short of fitting the data of Figure K-l(b). In an attempt to fit 
these data the following Hyper -Bet a Density Function was developed. 

Hyper-Beta Density Function . - Define a density function, similar to the 
Beta density as follows: 

P(x) • D x nb (1 - x 11 ) C for 0 < X < 1 (K ' 11) 

To establish D, it is known that the probability of a single event occurring 
somewhere in the region of interest is unity. Thus: 

/ p(x) dx s 1 

The constant D is found to be, 

nr (b+c+l+i-x 

D « — r— ^ 

r (b+£) r (c+i) 

This density, which has lower and upper limits of zero and one respectively can 
be transformed as the Beta Density was transformed. 
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Such transformation leads to the following density: 


„«.) - Da 1 - nb j4 ♦ ^ (1 
(00 , 


where : 


nr ( b+c+l+— ) 
n 

r(b+i) r(c+i) 



r(b+c+i+i) r(b+J) 

n ii- 

r(h+c+i+^) r(t+~) 

- r(b+H r(e+l) 

u n 

n r(b+c+i+^) 

n 


Comparison of this density for a set of constants b, c and n is shown in Figure 
K-U, again compared with the test data. Much better agreement is obtained 
because of the increased flexibility resulting from the* additional constant, n, 
in the definition of the density function. 
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Comparison of this density for a set of constants b, c and n is shown in Figure 
K-U, again compared with the test data. Much better agreement is obtained because 
of the increased flexibility resulting from the additional constant, n, in the 
definition of the density function. 



figure K-4. Comparison of the Gaussian, Beta and Hyper-Beta 

Density Functions with Test Data of Reference 12. 
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APPENDIX L 


Velocity Correlations . 


The auto and cross-correlation functions of the velocity terms can also 
be computed by use of the statistical flow model. These correlations are^ 
commonly called the Reynold’s stresses , and uv. Such ’’stresses 

produce a static pressure gradient which can also be obtained from the model. 
These correlations are developed below by application of the statistical 
techniques developed in Appendix J to the flow model. 

Axial Velocity Correlation . - The mean square axial velocity correlation 
can be obtained from the vortex flow field model and the statistical model of 
turbulence by the method defined by Equation L-l. 


“ 2 ■ N /n// 6m /a/ u(t)u(t) P (a) P (v 0m ) P(Y) p (n) dtdadv 0m dYdn 

(L-l) 


Where the axial component of fluctuating velocity, u(t), is Equation L-2: 

(L-2) 
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The respective vortex probability density functions are defined by Equations 
j-l4, J-15, J-16, and J-17. Substitution of the probability density func- 
tions and the expression for u(t) into Equation L-l and division by U“ yields 
the following integral: 
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Integration of this expression produces the desired mean square value of axial 
velocity fluctuation. 
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The last integral must be solved by numerical integration techniques. The 
vertical velocity correlation can be solved in a similar manner. 


Vertical Velocity Correlation . - The vertical component of fluctuating 
velocity, v(t), is given by Equation L-l and with the substitution X= U 0 t 
becomes: 


v (t ) = -nv 
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The mean square velocity will thus be: 
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Integration of this expression as far as possible in closed form yields the 
following equation: 
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Velocity Cross-correlation , 
components yields: 
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The cross correlation is computed, from the following integral equation: 
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Integrating this expression as far as possible in closed form yields Equation 
L-10. 
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Mean Velocity . - The mean value of the perturbation in axial velocity is 
found by use of Equation L-ll. 
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Substitution of the axial velocity perturbation and density functions into 
Equation L-ll yields the time averaged or mean velocity: 
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This can be simplified to Equation L-13* 
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The mean vertical velocity is found by similar methods. 
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After introducing the expression for the vertical velocity, this becomes: 
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And simplifies to: 


(L-15) 
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Static Pressure . - A static pressure gradient can be supported by the 
Reynolds shear stresses. This can also be obtained by use of the turbulent 
flow model . 
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Substitution of the density functions and integration yields: 
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The preceding velocity correlations and static pressure were evaluated 
by a numerical integration procedure for the probe location variation from 
h=o to h=H/2. The data are presented in Figure L-l, normalized by . 
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3.0 


Figure L-l 



. Velocity Correlations (Reynolds Stresses) and Static Pressure 
Difference as Computed from the Turbulent Flow Model. 
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Evaluations . - With the assumption of the vortex rotational direction 
being equally probable Plus or minus, and for the sensing probe in the center 
of the channel (h = H/2). The mean velocities and correlations can be found 
in closed form. These are: 
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At h = 0 (at one wall) the static pressure simplifies to: 
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A relationship between the velocity correlations (or Reynolds stresses 
as they are commonly called) and measurements taken by high response total 
pressure instrumentation can now be established by application of the turbu- 
lent flow model. An example is given in the next section. 
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Unsteady Total Pressure and Velocity Correlation . - The relationship 
between the turbulent fluctuation in total pressure can now be related to the 
fluctuations in velocity by use of the vortex flow model. This will link 
total pressure and hot wire anermometer measurements . The mean square velocity 
fluctuation was computed at the center of the duct in the preceding section. 
This is: 
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The mean square level of the total pressure fluctuations is given for 
the center of the duct by Equation L-19* 
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The ratio of velocity fluctuation to total pressure fluctuation is 
established from Equations L-l8 and L-19. 
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After appropriate simplification and for the velocity probability density func- 
tion having the exponents my = 4 and n v = lU (Refer to Figure 26 for the 
density functions), the ratio becomes: 
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Results of Equation L20 , which has very little dependence on the expo- 
nents of the density functions (mv> n v ) > were shown graphically in Figure 31 
d»f the main text for various levels of turbulence ( AP TRM s) and flow Mach 
numbers . 

This result is significant. For the first time a relationship has been 
developed between turbulence aB velocity fluctuations and turbulence as total 
pressure fluctuations. Previously, the relationship was obtained by assuming 
either sonic waves or a quasi-study analysis with a constant static pressure. 
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APPENDIX M 



FLUID DYNAMIC MODEL OF TURBULENT FLOW - COMFUTERROUTINE 

A computer program was written to evaluate the statistical properties of 
turbulent flow by numerical integration of the equations described in the main 
text. These properties include the Power Spectral Density function, Auto- 
correlation function, and Root Mean Square value of the total pressure fluc- 
tuations. Also included are the velocity correlation terms, u^, v 2 , uv, and the 

static pressure deviation, APg. A description of this computer program is the 
subject of this Appendix. 

The program is a digital computer solution of the Fluid Dynamic Model of 
Turbulent Inlet Flow. The single vortex total pressure variations is combined 
with the vortex random properties of size, strength, location, and spin direc- 
tion. The resultant equation is integrated with respect to the various random 
parameters and specified delay time (t) to yield the autocorrelation function* 
This autocorrelation based on the deviation of the total pressure fluctuations 
from the mean is made non-dimensional by the uniform stream dynamic pressure 
(qo = £ p Uo 2 ). The autocorrelation function is also computed normalized by the 
value at a delay ( r ) of zero (the mean square level). 

The Fourier transform of the normalized autocorrelation function is 
obtained by a numerical integration procedure. The result is the power spectral 
density (PSD) function of the total pressure fluctuations. 

INPUT: 

The input data card arrangement is shown in Figure 1. The input parameters 
are described below followed by a discussion of the Input default options. 

PARAMETE RS 

UO 
RHO 
HU 
H 
DP 

MO 
DTAU 

TAUL 


- duct (engine face) flow uniform velocity - ft/sec 

- duct flow density - LBM/ft^ 

- duct height (diameter) - in. 

- distance of sensing probe from lower duct wall - in. 

- ratio of root mean square total pressure fluctuation to average 
total pressure - AP^/^t 

- Mach Number 

- delay time increment for confuting autocorrelation function, 
normalized, DTAU = TAU * Uo/a 

- limit value for computing autocorrelation function 
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MVD, NVD 
MAD, NAD 
ANP, ANM 

LINE 

PN 

FLM 


vortex strength Beta probability density function (PDF) exponents 
vortex size Beta PDF exponents 

fraction of total vortices which have positive and negative, 
respectively, spin directions. Positive is counterclockwise; 
negative is clockwise rotation. 

Key for printing the vortex core size and strength PDF's. 

Number of proportional parts per decade of the frequency that 
PSD is computed. , /t _ T 

At- 10 1 '™ 

Limiting value of frequency that PSD is computed. 

1 g f g FLM -HZ 


INPUT PARAMETER DEFAULT OPTIONS 

IF Uo = 0, execution stops, use blank card to terminate. 

IF Mo = 0, Mo defaults to 0.4 

IF DP = 0, DP defaults to 0.02 

IF TAUL = 0, TAUL defaults to 3.0 

IF ANP = ANM = 0, set ANP = 0.5 & ANM = 0.5 

IF MVD = NVD = 0 and MAD = MAD = 0, then previous core size and strength PDF's 
are used, but not printed. 

IF MVD ^ NVD and MAD ^ NAD, then new core size and strength, PDF's are 
computed. 


IF LINE = 0, the PDF's are not printed, 

LINE / 0, the PDF's are printed. 

IF DTAU = 0, autocorrelation and PSD will not be computed. 
OUTPUT DESCRIPTION: 


The output data are printed in four groups - Probability Density Functions, 
velocity correlations, autocorrelation function, and the power spectral density 
function. Each output data group will be discussed separately. In addition, 
the input data are printed with each output data group. 


PROBABILITY DENSITY FUNCTIONS 

KV - Beta PDF constant for vortex strength (velocity) 

KA - Beta PDF constant for vortex size 

VBAR - Mean vortex strength (v 0 m fty /U o ) 

ABAR - Mean vortex core size (a/H) 

A/H - ratio of vortex core size to inlet duct height 
A - vortex core size - in. 
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P(A) - probability of vortex having size A/H 

VTM/Uo - ratio of vortex maximum velocity to local flow velocity 

(vSnax/Uo) 

VTM - probability of vortices having strength VTM 


VELOCITY CORRELATIONS 


H/HU - 
H 

UUBAR - 

WEAR - 

UVBAR - 
UBAR/Uo 
VBAR/Uo 
DP/Qo - 

DUUDY - 

DWDY - 
DUVDY - 
DUDY - 
DVDY - 
DFBY - 


total pressure probe location; h/H in model 
absolute value of probe location - in. 

mean square axial velocity fluctuation, normalized by local 
velocity square (Uo^) and the term (NH/Uo) . 

mean square lateral velocity fluctuation, normalized by Uo^ and 
(NH/Uo) 

mean velocity cross correlation, normalized by Uo^ and (NH/Uo). 

- mean axial velocity fluctuation, normalized by Uo and (NH/Uo). 

mean lateral velocity fluctuation, normalized by Uo and (NH/Uo). 

mean static pressure fluctuation, normalized by dynamic pressure 
based on local flow (qo » •§■ p Uo*) and(NH/Uo). 

gradient of mean square axial velocity fluctuation with respect 
to Y. 

similar to DUUDY for mean square of lateral velocity 
similar to DUUDY for velocity cross correlations 
similar to DUUDY for mean axial velocity 
similar to DUUDY for mean lateral velocity 
similar to DUUDY for mean static pressure 


AUTOCORRELATION FUNCTIONS 


ALUMP - intermediate results 

ABAR/Uo - inverse of time delay normalizing factor 
AK - number of vortices per unit time (N) 

TAU*U0/A - normalized time delay 
TAU - actual time delay, t - sec 

RXT - actual total pressure autocorrelation at time delay, r 
RXT/RXTO - total pressure autocorrelation normalized by value at r= 0 
RUT - actual velocity autocorrelation at time delay, r 
RUT/RUTO - velocity autocorrelation normalized by value at r = 0. 
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POWER SPECTRAL DENSITY FUNCTION 


F 

FR 

GXF 

GXFR 

RMS 

GUF 

GUFR 


frequency HZ (CPS) 

normalized frequency FR = F* ABAR/Uo 

total pressure power spectral density at frequency, f 

normalized power spectral density at frequency, f orFR 

integrated area under PSD curve 

velocity power spectral dens it y at frequency, f 

normalized power spectral density at frequency, f or FR 
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FIGURE Ml INPUT DATA DECK ARRANGEMENT. 
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PROGRAM 

Ok 815 TRAOF 

00C 6610 FTN V4.0-P332 OPT=l 

03/15/73 10 .12,44. 

PAGE 

1 


PROGRAM 0 481 5( INPUT ,PJTPU T , T fl n E5 = IN 3 UT ,T AP E6 =0 U TP1J T ) 



C AUTOCOP r l AT I ON PROGRAM 

VEPSION TWO BEGUM 11/30/1972 

A 

15 





RR AL INTVA, INTVR, 

INTV0 , ’’NTT A, TNT T9* TNT TO, I NT A, INTB, 

A 

20 



ro — 

o 


lINTC,INTn,TNAA,INAP, 

INA 0 ,KV,KA 

A 

30 



*= — 

5 

R5 Al_ INTYAjINTYR, 

INTYO 






DIMENSION TR(IQO), RXTT(IOO), FUN ( 20 ) ,RUTT(100) 



COMMON 

APAR, 

AK, 

ALUMPA, 

0 A , 

0T, 

DT AU, 



HU, 

A 

60 


1 

INTVA, 

INTVO, 

KA, 

<v. 

LIME, 

M, 

MAD, 

MV, 

A 

70 


7 

MVO, 

Nj_ 

NA , 

NAO, 

MV, 

N VO* 

RHO, 

T_L» 

A 

80 

ir 

3 

uo, 

VRAR, 

MA , 

TAUL, 

ANP, 

ANM, ALUMP, 

1 'FUN#>: ; 

A 

90 


4 INTVH 

COMMON /OAT Af TR,RXTT,PN,FLM,TLM,NT,RUTT A 100 

10 PQpHAT <1H1,2QX,23HAUT0CQPREL "TION PROGRAM//) A 110 


2 0 5QPHAT (8F10.5) a 120 


15 

3 0 FORMAT ( 14X , 2HU0 «5X*3HRH0*4X * ?HHU* 5* * 1HH , 5 X3 HQ P T* 4 X * 2HMO , 3X , 4HQT A.U__ 

A 

JULG 


1*2X*4HTAUL/10X,F7.1*F 7 .4*2F7 . 2,F6 . 3, F7 . 3 , F6. 3 , F7. 3 > 

A 

140 


4 0 FQPHAT (15X*F7. 3.F8.4.F11 ,5, F c ,5.F11.5,r9.5) 

A 

150 


FORMAT </15X,8HTAU*UP/A*3X,3HTAU.7X.3HRXT.4X,8H l? XT/RXT0, 

A 

160 


* 5X,3HPUT ,4X ,8HRUT/RUTO/1 

A 

165 

20 

60 FORMAT ( / 1?X ,2HMV, 3X ,2HN V ,6X , 2HK V * 5* , 4 HVRft R/ 9X , 215 , F9 . 0 , FR • 4//12X , 

A 

170 


12HMA, 3X,?MNA*6X,2HKA,5X, 4 HA R A p , 5 X , 2H N* ,5 X , 2H N- / 9X , 21 5 , F 9 . 0 , F8 . 4 , 

A 

180 


* 2F7.7) 

A 

185 


70 FORMAT (415,2510.6,215) 

A 

19t 


*0 FOPMAT </12X, 5HALUMP, 7y, 7HA3 AVU0,7X, 7H A< /10X,3E12.4) 

A 

200 

25 

C AWHOL c RUNPH OF CONSTANTS 

A 

210 


CA=?. 0 

A 

220 


CB=SQRT (2*0/3.01 

A 

230 


CC = SORT(2.nj ... ____ 

A 

2 40 


''9 = 5.0/16.0 

A 

250 

30 

CE=li.C/(16.1*CC) _ _ 

A 

260 


r F = SOBT <6 * 0 ) 

A 

2 70 


A LUMP A=EX°(1.0) /CA 

A 

280 


ALUMPB=FX°C1.5) *CP/9. * 

A 

290 


ALUMPD=EXP(2.0) /CC/4.1 

A 

3 00 

35 

90 R c AO (5,29) UP . RHP . HU, H , nT, T L , PT AU * T AUL 

A 

310 


IF (UO.EO.C.O) STOP 

A 

320 

IF (TL.EQ.0) TL=0.4 

TF (OT.EO.O.J) nr=c,fi 


IF (TAUL.EO.O) T AUL = 3, 0 

A 

340 

40 

OOTAU=2.0*OTAU __ 

A 

350 


P c AO (5.79 MVO* NVO. MAO, NAP, AMP.ANM.LTNE 

A 

360 


TF f ANP,N5, 0,0) GO TP 95 

A 

362 


IF (ANM.NE, C. r > GO TO QR 

A 

364 

AN P = 0 .5 ... __ . 

45 

AN V = 0 • 5 




3 5 R c en (5,20) PN,5LM 

A 

370 


IF (MAC.EO. o .ANP. NAP. £Q. -> GO TO ICO 

A 

380 


GO TP 110 ____ ... _ 

A 

390 


ICO 15 (MVD.5Q. O.ANO.MAn.EQ. " > GO TO 120 

A 

4C5 

30 

11 0 CALL FROFLE 

A 

4 1C 


120 L T *' E = 1 4 

A 

420 ... . _ 


C L r ADT NG COMSTAN SLUM* 

A 

430 _ ... ......... 


AL (IMP = 8 . 0 *KA*HU/(12 .G*UC) _ 

_.A 

44C_. ... 


CALL VELCOR ...... 

A 

.45. C . . ... ... _ . .. .... 

55 

IF LOTAH, rc. 0.?) 50 9^ . . _ 


, 


PROGRAM DU815 T P ft CE HOC 6600 FT N V4.0-P33? OPT=l 03/15/73 10 .12.4 4. PAGE 2 

ALUHUA^AK »KA»HU/(12. Q»UO) »EYP(1.?) ' ~~ 

TLP=flRAR*HU/<12,q*UQ) A 460 




TLLM=1. 0/(6. 203185 

*FLM) 

A 470 



tulm=tllm/tlh 


A 48C 

60 


CCC=10.0E-08 


A ■ 490 



DO 140 1=1,10 


A 500 



IF (CCC-TULM) 130, 

15(3, 15 0 

A 510 


170 

ccc=io.o*ccc 


A 520 


140 

continue 


A 530 

65 

150 

DTU=0 • 15CCC 


A 5 40 



IF (DTU* GT • DTAU) OTU=DTAU 

A 550 



TU=0 • 0-DTU 


A 560 



WRITE (6,10). 


A 570 



WRITE (6,30) UO,RHO,HU,H,OT, tl,OTAU,TAUL 

A 580 

70 


WRITE (6,60) MV , NV 

,KV, VBAP,MA,NA,KA, APAR,ANP,ANM 

A 590 



WRITE (6,80) ALUMF 

, TLM , AK. 

A 600 



WRITE (6,50) 


A 610 



DO 240 ITU=1,1000 


A 620 



IF (TU.GT.OTAU) GO 

TO 160 

A 630 

76 


TU = TU +OTU 


A 6 40 



IF (ITU*EQ*1) GO TO 1? 0 

A 650 



DTU=2>0TDTU 


A 660 



GO to 170 


A 670 


160 

IF (TU.lT.OOTftU) TU=OIflU 

A 680 

80 


TU=TU+DTAU 


A 690 


17 0 

TR-(ITIJ) = Tij 


A 700 



IF (Til • GT * T AUL) GO 

TO 25 0 

A 710 



TAU=TU*TLM, 


A 720 



tau?=tau*tau 


A 7 30 

*6 


AOH = 0.0. , 


A 740 



J= 0 


A 750 



FTA^o.O, , . 


A“ 760 

F,I P = 0 . 0 . 



FIP=0,0 , 


A 780 

9 r 


FIT A= 0 1 0 . . ; 


A 790 

■ FTIB=%0 . .. 



FI TO= 0*0 


A 810 



DO 20 ) .1=1,2 r 


A 820 



ftOH = A DH + n A 


A 830 

96 


IF (AOH.GT.l.Q) GO 

to 21 * 

. A 840 



ft= A OH *HU ; 


A 850 



pah=-h/a 


A .8 60 



RAH2=RAH*RAH 


A 870 



RAHU= (HU-.H) /ft 


A 88C 

g IOC 


RAHU2=RftHU*FAHU 


A 69C 

U1 


PH I AL =-R AH 


• A 900 



PH I AH = RAH U 


• A ; 9 If 



KOnr=l 


A 920 



CALL SIIHMEP (PHIAL 

,CIM"TAL,W0PE) 

A 930 

106 


CALL SUMMER .(PH I AH 

,CIMTflH,KOP r ) 

A 94C 



TNTYA = <RAH*FXP< -RAH2)-RAHU*E XP (-R A H J 2) + CINTAL+C INT AH) 

A 950 



INTY9 = EXP <-1.5*RAHU2)Mi . 0-3 .0*RAHU2)- r XF(-l ,5*PAH?) * 

(1*0-3. *RAH2> A 960 



CHIL= -CC*R am 


A 970 



CHIHsCC^RAHU 


A 980 

11? 


CALL SUMMER (CHIL, 

CI^T DL , KOD E> . , 

A 990 
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CALL SUMMER (CHlH,CTNTOH,KOO r > 

A1000 


INTYO^EXP (-2.0*PAH2> *( RAH*R A H2 /4 • O+CO* RAH) -EXP (-2. 0*PAHU2> * (RAHU*R 

A1010 


1AHU2/4 • 0 + CO*RAHU> +CF* (CINTOL+CINTOH) 

A1Q20 


FUNnA=FUN(T) 

A1030 

JY> 115 

UOA=12.0*'JO/A 

A 1 0 4 0 

— r~ ■ ■ 

U0A2=U0A*U0A 

A 1 J 5 0 


TAUA-=TAU*TAl)*UOA2/4.0 

A1060 


TAUR-TAUA*4. 0/3.0 

A 1 0 7 0 


TAUD-TAU2*UOA2/?.' 

A1080 

120 

TA = C X°(-T AU A ) 

A 1 0 90 


TP = EX° C-T AUP> 

All 00 


TO=EXP C-T AUP) 

Al 110 


Tt=0.0 

All 20 


PHI=UOA*TAU/2.0 

All 30 

125 

PSI=UOA*TAU*CB 

_ AU40 


0M€GA=UOA*TAU/CF 

A1150 


CHIrUOA*T AU/CC 

A 1 1 60 


KODE= 2 

Al 170 


CALL SUMMER (°HT , TNT A, KOOE) 

Al 1 80 

13 P 

CALL SUMMER fPSI • TNTP, KOOF) 

Al 190 


CALL SUMMER (OM r G A . INT C , KOD* ) 

A1200 


CALL SUMMER fCHI , INTO, KOOH 

A1210 


IMTTA=INTA*Tfl 

Al 220 


INTTB=TP» (INTR+TNTC) 

A1230 

135 

INTTD=INTO*TO 

Al 240 


P A A = FIJNCA *If!TTA* INTYA 

A1250 _ _ . 


PAP=FUNCA*INTTB»INTVp 

A126C ..... __ ....... 


pao=func A ♦TNTTD^ lNTYp 

A127C 


J=J*1 

„ A 1 2 8 0 _ 

140 

GO TO (180,190* , J 

Al 290 


lf»n FI A = FIA-*-P AA 

A13GC 


FI0=FIP+PAP 

A1310 


FIO-FID+PAO 

A1320 


GO TO 200 

A133G 

14* 

1°0 J= 0 

A1340 


FIIA=FHA+PAA 

Al 350 ____ 


FT TR=FIIB+P AR 

A 1 360 


fiio=fitoffac 

A1370 


?rn CONTINUE 

A13BC _.. ... _ _ 


210 CONTINUE 

Al 390 


IN A A = (FAA>4.n*FTA*2.C*FTTA) * PA/3. 0 

A140C ..... _ 


TMAB = (PA 9+4 .0*FI°+?,3*FTIB) ? r A/3.C 

A141C 

TNA9*TN~AP*< ANP-ANM) 


INAO=(FA04-4.0*FTP+?.r*FlTP) * PA/3.0 

A1420 

155 

<?X T = A LUMP * ( ALUMP A*IN A A *INTVA+ UU MOP* IN AS *1 NT VP ♦ AL ! J V °P* I N AO* IN TVO ) 

Al 4 30 _ 

RXT = RXT*A1< 

®UT=flLUHUft»TNAfl*IMTVA 


TF fITU.NE.il GO * r O 2^0 

A1450 


pxjn=Rx T 

A1460 

16' 

R.U.T_0=RUT 



220 PXPTrRXT/RXTO _ 

Al 4 7 (L 

®urt=put/ru t o ... . 


PXTTf ITU) =RX°T 

A 1480 ... .... 

PUTT ( ITU ) -°IIPT 

165 

L INE=LIN C «-i 

Al 490 



— - - - 

- -• - - - — 

- 1 




pro gram 


0^815 


TRAC^ 


non 5 600 FT N V4 « 0-P3 32 0PT=1 03/15/73 10,12.4 4 



IF ( L INE . LT , 61 ) GO TO ?30 

A 1 5 0 0 


WRITE 16,10) 

A1510 


WRITE (6,50) 

A1520 


L I NE- 6 

A 1 530 


2T0 HR ITT (6,40) TU , TAU , FX T , RXRT ,RUT , RURT 

A154C 


IF (PXRT.LT. 0.001) GO TO 250 

A 155 0 


2un CONTINUE 

A1560 


2E 0 CONTINUE 

A1570 


NT = IT U-l 

A1580 

5 

CALL POWER 

A1590 


GO TO 90 

A1600 


r ND 

A161C- 


SUBROUTINE FA C T 


TRACE 


CPC 6600 FT N VU.0-P332 OPT = l 03/15/73 10 .12.53 


SUBROUTINE FACT (TS.SU M ) 8 10 

SUN-1.0 B 2C 

DO 10 1-1 .IS 8 30 

SUM = I»SUH B 40 

10 CONTINU E B 50 

RETURN B 60 

B 70- 


SUBROUTINE KXB&9 


TRACE 


SUPROUTINE KXQA R (H,N,K,XPAR) 

REAL K 

IS=N+N»1 

CALL FACT <M«SKM ) 

CALL PACT (IS.SKNN) 

CALL FACT, <N,SKN) _ 

KsSKNN/<SKW»SKN) 

IS=IS+1 

SKNN=5KNN*I5 

.. SK*=*KH»(H+n 

y3AR=K»SKH*SKN/SKNN ’ 

RETURN 

ENP , 




SUBROUTINE POWER 


me? 


CPC 5600 FTN V4.0-P332 QPT=1 03/15/73 10,13. 32 


PIGE 




SUPROUTTNF POW c R 

0 10 


DIMENSION TRI100). RXTT(IOO) . EE(IOO). EU(iOO), RUTTflQO) 

COMMON /DATA/ TR.pXTT.PN.FLM .TLM.NT.RUTT 

IX) 


10 FORMAT (1H1.15X.19MF ONER SFECT«>UM HACK /I 

D 40 


s - 

c 

20 FOPMAT (1H0.16X.1HF.OX.2HFP. 9 ''.3HGXF.BX. AHGXFP . flX, 3HPMS, 

D 40 




* 9X. 3HGUF. 9X. 4HGUFX/1 

0 45 




70 FORMAT (10X.FF12.4) 

0 60 




40 FORMAT (F12X.5H PN , FX, 7HAB AP/U0.7X *FH FL IH IT/10X . 3E12 .4) 

O J.Q 




F=1.0/PN 

0 80 



If 

AHUL=10.0**F 

0 90 




WRJTF (6.10) 

0 100 




IF (FlM.EQ.O.n FLM=1 3 • G c +f 6 

0 lie 





0 120 





D 130 



l c 


0 140 




F=1.0 

D 150 




EEfl>=0.0 

O 160 


EU ( 11 =0 • 0 



DO 50 1=2. NT 

O 170 



20 

EE (I) =ALOG(PXTT(I) ) /TR{I> 

O 180 


rum = ALOG(PUTT(in/TR(I) 



50 CONTINUE 

O 190 




RMS-0 ■ 0 _ .... _ 

O 2QQ 





O 210 . ... ... 



?5 

IF (1-2) 80 , 6 P • 7 0 

0 220 ... . 




63 F= AMUL 

.. .... D 210 . 




50 TO 80 

D 240 




7Q F=F*AMUL 

D 250 




*0 CONTINUE 

0 260 



*n 

IF (F.GT.FLM) GO TO 150 

0 270 




fp=f»tlm 

0 280 . .... _ 




ARG=2.0*3.i4159*FR 

0 290 




TRPONS=l. OAF R 

0 300 




TF (TPCONS.GT.TP (NT) ) GO TO «3 

O 310 



35 

OT = 0. l*Tp0ONS 

D 320 




GO TO 130 

0 330 




P-0 OT=3.1*TP(NT) 

0 340 




ICO AN = TP(NT)/DT 

O 350 




N= AN<M _ 

D 360 . 



Jtl 

TA = TF(l)-OT 

O 370 




GXF=0 .0 

. . D 360 . . . 


GUF=fl . 0 . . • 



00 123 J=1.N 

O 390 




T A=TA +OT 

n 400 



4 c 

CALL TPLU°1 (TR. ctr .TA. ry . PFXOTP.NT) 

D 410 


rftlL TPLUP1 -fTP.cn.Tf . r7.nrftn T o.MT) 



RXT=FXPf FX*T A) ...... _ . 

fl. 420 ..... 


PllT=FXP(rq#TAl 



IF (J.FO.l) GO T 0 110 _ 

0 430 



Rf 

T - T A— 3 . 5*OT 

n 440 




ANF-C OS ( A RG* T) 

D 45.0 _ 




GXF=GXF+(RXT4.RX t L)*ANF»OT*2. ■: 

D 460 


GUF=GDF> 1 PUT ♦RU T LJ *4Nf *DX*2. 0 .. _ ...... 



_J.in pxt.L-RXT _ . ... ... .. .. _ 

0 470 



55 

RUTL=RUT 




_ _ _ __ j 


'1 



SUBROUTINE 

pnwc 

:R TRACE 

COC 5 6: H FT N V4.C-F332 OFT = l 

13/15/73 1C. 13. 32 



1?0 

CONTINUE 


0 

400 



IF (I.EQ.l) GXFF=GXF 


0 

49C 

IF (I.EO.l) GUFF=GUF 



GXFR-GXF/GXFF 


D 

500 

SUFR=GUF/GUFF 



IF (I.EO.l) 'GO TO 13 C 


0 

510 



RMS=RHS+0 .5*(GXF+GXFL> 

MFR-FRL) 

0 

520 


130 

COM T INUE 


D 

530 



WRI T E (6*30) F,FR,GXF* 

GX FR* RNS *GUF * GUFR 

D 

540 



GXFL=GXF 


0 

550 



FRL=FR 


0 

560 


14Q 

CONTINUE 


0 

570 


150 

CONTINUE 


0 

500 



RETURN 


D 

,590 



ENT 


0 

600- 


SUP POUT TNF PEOPLE 


TRACF 


PC 6600 FTN V4 . Q-P332 OPT^l 03/15 /73 1(1.13*09 


PAGE 


I 




SUBROUTINE FROFLF 

E 

10 



DIMENSION FUN (201 

E 

20 

-fQ- 


REAL TNTVA,INTvn,KA,KV,INTVB 

E 

30 

H 

tO- 


COMMON ARAR, AK,ALUMPA, DA, DT, DT Ally H, HU, 

r 

4r 


5 

1 INTVA, TN T VD , KA, KV, LINE, M, MAD, MV, 

E 

50 



2 M VO, N , NA, NAD, MV, NVO, PHO , TL, 

E 

60 



* uo, vpar, ma, taul, anf, anm, alump, fun. 

F 

7 C 

4 I NT VP 



10 FORMAT (1H1, 2CX,25HPP0BARTLI TV DISTRIBUTIONS/) 

£ 

80 


If 

20 FORMAT (14X , 2HU0,5X, 3MRHO ,4X , "HHU , 5X , 1 HH , 5 X3 HDPT, 4 X , ?HMO , 3X ,4HDT AU 

E 

90 



l,2X,4HTAUL/lt)X,F7.1,F^.4,2F7. 2,F6. 3, F7 .3 ,F6. 3,F7.3> 

E 

100 



30 FORMAT (/12X,2HMV,3X,?HNV,6X, 'HKV,5X, 4HVBAR/9X, 215, F9. 3, F6.4//12X, 

E 

lie 



12HMA, 3X,2HNA,6X,2HKA,5X, 4 HA 9 A? ,5 X , 2MN* , 5 X, 2HN- /9X , 2 15 , F9 . 0 , F8 . 4 , 

E 

120 

* 2F7 • 4) 


15 

4 a FORMAT </ 11X , 3H A/m ,5X, 1HA,5X ,4HP<A) , 3X ,6HVTM/U0 ,4X , 3HVTM , 4X ,6HF ( VT 

E 

130 



1M) ) 

E 

140 



F0 c OPMAT (8X,F6.?,F7.2,F8. 4,FR .? ,F8.2, FP .4) 

E 

150 

POA = 0 .« 



MV=MVD 

E 

160 


20 

NV=NV9 

E 

170 



M A = MA Q 

e 

180 



NA=NAO 

c 

190 



CALL KXBAR ( MV, N V, KV , V «A R ) 

E 

22C 



CALL KX9 A R (MA,NA,KA,ARAR) 

E 

230 


25 

NPPNT=tIN2 



IF (NFRNT.FO.O) GO TO 61 



WRITE (6,10) 

E 

200 



WRITE (6,20) UO,RHO,HU,H,OT, TL,DTAU,TAUL 

E 

210 



WRITE (6,30) MV,NV,KV, VBAP,MA,NA,KA, APAR,ANP ,ANM 

E 

240 


30 

L I ME= 1 1 

E 

380 



WRITE (6,4ft) 

E 

390 



FI M = M A ♦ ? 

E 

250 



IS-MV+2 

E 

260 

VA1 = ( T c -1 ) * IS 


35 

IS0=MV+NV+3 

c 

290 

vAP=(ISD-l)* xso 

intva =v a i /van 



V A 1 = V Al* ( IS+ 1) 

z 

320 



VAO=VAO»<ISO+1) 

F 

330 


4T 

IN TVR^VA 1 / V AO 





VA1 = VA1MIS+?) 

E 

350 



V Apry A D* ( IS 0+2 ) 

. E 

360 

INTVO = VAl / V A 0 " 



SV = KV 

w 

4 0 C 


45 

fa~=ka 

c 

410 



DR=0. 11 

F 

42ft 



AR=0. 3-OP 

E 

430 



AMAX= 3,n 

c 

440 



N= n 

sr 

45C 


6 0 

OO 103 T = 1 , 1 31 

E 

460 



AR= AR+OR 

z 

470 



A= AR’HU 

z 

480 



V r H=AR*UO 

E 

490 

1 


IF GAR.NE.O. 1) GO to 6ft 

r 

500 


6> 

I r (MA.EO.O) FOA=SA 

E 

510 
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£ PRQFLE TRACE - 




TF 

fMV.EO.P) FOV=SV 




GO 

TO 70 



60 

IF 

fAR.GT.l.O) 0 0 TO 110 




Fnv=SV*ap**MW*<1.9-AR> **NV 


• P0A=SA*AR**MAM1. O-AR) **Nft 


70 

AHAX=AMAX1( ABAX , PQA) 




IF 

CI.EO.l) GO TO 80 


RAT=POA/ANAX 



IF 

<R AT ,LT • 0 ,00 01) GO TO 80 




N= N4-1 



PO 

LINE=LINE+1 




IF 

(N°RNT.£Q.Q) GO TO 100 



IF (I INE« LT » 61) GO TO 90 


WRITE (6.10) 
WRITE <6,40) 


LINE = C > 



°0 

WRITF (6,50) AR, A,POA, AR,VTH,POV __ 


iro 

CONTINUE 


110 

CONTINUE 

5 


RETURN 
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SUPROUTINE SUMMER <X,£RF,KOnP) 

F 10 





r 

NO. 70804 IN THE H-P PROGRAM ROOKLFT 






H 

DIMENSION XT C 30) 







DATA (J=l) 

F 30 





5 

50 TO (10,30), J 

F 40 






10 DO 20 K=1 , 3 G 







N= 2* C K-l ) +1 

F 60 






E=303,0/N 

F 70 






20 XT(K)=13.1**E 

F ftC 





10 

J= 2 

F 90 






CHFCK=0. C 3 0 C 1 

F 100 






30 SUM^Q.O 

F 110 




TF (X.LT.2.5) GO TO 40 



ERF=1 .0 

F 130 





15 

GO TO 70 

F 140 






lp D=1,0 

F 150 






A N = 1 • 0 

F 160 




DO 50 1=1,30 



IF (X.GT.XT (7)) GO TO 60 

F iflO 





2C 

N= I-i 

F 190 






AN = 2. ft*AN 

F 20 0 






IF (N.EO.0) AN= 1 * 0 

F 210 






NO=2*N+i 

F 220 






D=P*ND 

F 230 





25 

XX=X**MP 

F 240 






T EPM= AN*XX/P 

F 25P 






IF (TERH.LT. CHECK) GO TO 60 

F 260 






SUN=SUM+TERM 

F 27C 




R=TERH/SUM 


30 

IF CR.LT. CHECK) GO TO 60 







5" CONTINUE 

F 280 






F0 continue 

F 290 






FRF=1 . 1263ft* EX P(-X*X)*SUN 

F 300 






GO TO (00,60), KOPE 

F 310 





35 

PO ERF=j ,0-ERF 

F 320 






°0 ERF=ERF/i. 1283ft 

F 33C 






RE TURN 

F 340 





■£HT_ - F 350- 
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COC_6 600 FT N V4.Q-P332 QPT=1 0 3/ 15/73 10*13.30, 


SUFPQUTTNE V ELTOR 

DIMENSION FUNC2Q) ■ 

DIMENSION DUC25) ,nUU<?5) ,QV( ?5) ,0V VJ_25 ) , OUV( 2 5 ) , QF F t 25) 

REAL INTV A ,INTyn,KA,XV,INTV B 

COMMON ABAP, AK,ALUMPA, ' OA , OT, OTAU, H, 

1 INTVA , IN T VO , KA, KV, LINE, M, MAD, 

l MVP, N> NA , NAD, MV, NVO, PHO, 

AMP, ANH. 


.■? UP, VPAR, MA, TAUL, AMP, ANM, ALUHP, 

4 I NT VP 

30 FORMAT (1H1,2QX,23H VELOCITY CORRELATIONS//) 

70 FORMAT ( 14X , 2HUO ,5 X , ?M R H0,4X ,?HHU, 5X ,1 HH ,5X3 HDPT,4X , 2HM0,3 X 

1>?X»4HTAUL/19 X, F7 . 1 , F7 . 4 . 2F7 .2.F6.3.F7 .3.F6.3..F7.31 ~ 

POPHAT C/l?X t 2HMV,3X,?H NV,eX , ?HKV,5X ,4HVBAR/ 9X , 215 , FR. 1. F9 . 

12MMilXf_2HNA,6X,?HKA,5X,4HABPq,5X,?HN+,5X.2HN-/qX,2I5,Fq.O, 

* ?E7 . 4) ~ — — 

U* FO RMAT (5X,F10.? , F8 . 3, FR « 5, 5 F * ,g) ‘ 

c O FQPHAT ( / 11X ,4HH/HU ,5< ,1HN, 5 X .5HUUPAR, 4X ,5HV VB AR,4 X , 5HUV8AR 
19 AP/UQ, ?X , 8HVRAR/UQ ,1X,7H ne/QO/) 

EQ EQPHAT <1H0,15X,?QH( U/ UP ) 2/ ( P r T/QQ) 2 * ,F7.4 ) 

65 FQPHAT (/ /I IX, 4HH/HU ,5 X , 1 H H , 6X, 5H0U UbY ,4 X , 5HOV VOY, 4 X t 5H0UVD 
» EX, 4H0UDY ,5X,4H0VQ Y. 5X , 4HD' 3 CY>) 

EQ UIVALENCE (DPTPT6,PT). (AMO ,TL) “ ~ 

AHAsMA . 

AN AaN A ’ _ 

AN V=H V ' 

ANVsNV 

E AT = K4 * V9APMAMV + Z. n )/ ( A M V-f S N V+3, 0) 

ALUMUUs?, 405»FA T 
A L II MV V = ? , 4 0 5 *F A T 

ALUMUV=-1.7sq»FAT 

A LUNUs4.13?74»KA»VPAP» (ANP-A NM) 

ALNMVsB. 2R744 *KA*vsAP» ( A N P- A N > 

ALUNQ P= 4.81 BQ3»FAT " 

AN = N * 

Q-A=AN/20Q 0 «_Q _ “ ~~ 

UUCPM = 1. 3/(4, 0±.i ,3 . 6eS2* < AlLVijuiL) »XAlMV.±3_ ^_) / C ( AMV+ A NV+5 . 0) * 

DLt.Cu.0_)J_> 

Q°ia=_l.Z? AM&*AMOZJ (_ i_,j 4-JL . 2 * A m r* ft M O l _*Z3_,Ji ___ 

ANNsC OPT PTP/Q PTR) »»2 

OU O=UlJO°TO » ANN I. .Z I 

ANHUO -U'iOZ?. , .a* ( A N y.t a [O v + ? , C L * ( AN V > ANV +3, Q \ 

IXAil A.± A N_A ♦ 2^_fl 1 * X At! A + AJLA.+. !.■ . )/ ( 8.53Q73»( AMV» 1,0) * ( *HV + ?,fn 

_* .. -M.AMA±:UJUZ< A.NA + 20)) 

AK = L.0_.. ... . 

H-RIlE-.J64.iOJL . . .. . „ " 

JLRITE--1A4 20J__U0.,.RH0, 6U , H ^O.T.,.T_L.,0 TAU, t ft qlI ‘ 

H_RItE_L6,3PJ. MV3NV,ZV4.mA.P^MA.,m4KA4APAR^^{P_,AtiN ’ ‘ * 

-HR1!£ .. 16^51 1 

.. DH.P_=_*_Q_5_. . _ ... 

RS.-_Q « 1-J3H® 

00.12 ZJ = 1, ? onn ' _* ^ " " 

HRrMR+PHa 

IE. t_MR,.G.L,_C. 5, 1L r ? o„ jn 171 

- - IZ=I„ " " 

HI=- HR» HU _________ * 


4//12X, 

Ffl.4, 


H 210 
_^hL_2ZC_ 

H 9_91_ 

— H 9.0H. 


H 2 fer 
_-iL2>(L_ 
._ 1 ±. 280 . . 
—H.Z9.C 


SUBROUTINE V r L COR 


TRACE 


CDC 6600 FT N V4.0-F33? QPT=1 03/15/73 10.13.30 


AOH^O.O 

H 

340 


J=0 

H 

35C 


FIU=0 • 0 

FIV=0.0 

FIP=0.P 

F I UU= >3.0 

H 

360 


FIVV^O.O 

H 

370 


FI UV= 0.3 

H 

380 



FI IU= 3.0 


ftiv=3.o 


FI IP- 0 • 0 


FI IUU=0 • 0 

H 

390 


Fiivv=o.r 

H 

400 


FIIUV=0.0 

H 

410 


^ 00 10? K= 1,2000 

H 

420 


AOH=AOH*0A 

H 

430 


IF f AOH.GT *1.0) GO TO 110 

H 

44C 


A= AOH*HU 

H 

450 


R AH--HI/A 

H 

460 


; RAH2=RAH*R AH 

H 

47P 


RAHU= (HU-HI) /A 

H 

480 


RAHU2.=RAHU*R AHU 

H 

490 


PHIAL=-PAH 

H 

500 


PHIAH=RAHU 

H 

510 



SI = 0. 7Q7*(HU-HI) /A 


SII=Q .7D7*HI/A 


KOPE=l 

H 

520 

CALL SUMMER (PHIAL,CI'ITAL ,KO nr ) 

H 

530 

CALL SUMMER fFHI AH,CIM TAH ,KO 0 = > 

H 

540 

i CALL SUMMER (SI , SIN, KODE 1 

CALL SUMMER (SIT , SIIN, KOOH 

IF ( I .GT • 11 GO TO 70 

H 

550 

F2=A0H**M 

H 

560 

F3= ( 1 * O-AOHI **NA 

H 

570 

! FUNCA=F2*F3 

H 

580 

FUN (K > =FUNC A 

H 

590 

70 FUNCA =FUN ( K 1 

H 

600 


PAU=FUNCA*(FXP(-.5*RAH2) - EX 3 (-“.5*RAHU2> > 
PAV=FUNCA*(SIN+SIIN1 ~ 



PA UU= FUNCA* ( RAH * EXP (-RAH21-R AHU* EXP ( -R AHU2 ) + CINTAL + CINT AH) 

H 610 



PAVV=FUNCA* (CINTAL + CI'ITAH) 

H 620 



PAUV= FUNCA* (EXP (-RAH2) -E X P< - PAHU 21 1 

H 630 


PAP=P AVV 


J = J + 1 

H 640 



GO TO (80,901, J 

H 650 



PO P I UU= PAUU + FI UU 

H 660 



FI VV=PAVV*FIVV 

H 670 



FT! IV = P AUV ♦ F I UV 

H 680 


FTU=FTU+PAU _ 


FI V=FI V*P AV 



FIF = F IF+F AF 


GO TO 100 

H 690 



°0 FI IUU = PAULI* FHUU 

H 700 



FIIVV=FAVV+FITVV 

H 710 



FITUV=PAUV*F IIUV " H 720 



T RACE 


SUBROUTINE V r LCOR 


CPC 06 *" FT N V4.T-F332 OFT=l 03/1 5/ 7 7 IP, 13. 30. F AGE _3 


FIIU=FIIU+PAU 

FIIV=PIIV+PAV 

CU- 




FIIP=FIIP+PAP 




M 

a> 




J= 0 

H 730 




115 


ICO 

CONTINUE 

H 740 






110 

CONTINUE 

H 750 







S I NUU= CF AUU* 4. 0 *FIUU+2 . a * FT T UU ) * DA/3 . 0 

H 760 







SINVV=(PAVV+4. 0*FIVV+2. 0*FII VV)*OA/3.0 

H 770 







SINUV=CPAUV* 4.0 *FIUV+2.n*FIIU7)*0 A/ 3.0 

H 78C 




120 



SINV= (PAV + 4. 0*FIV+2.0*FIIV> *HA/3.0 




SINU= (PAU+4 . 0*FIU+2. C*FI Ilj) * D 5/3 . 0 

SINP=CPAP+4. 0*FIF+2.‘ »riIF> *01/3. * 

ubarr=ai_umu*sinu*ak 

vbapr=alumv*sinv*ak 


125 



DP=ALUMOP*SINP*AK 




UUFAR=ALUMUU*SINUU*AK 

VVPAR=ALUMVV*SINVV*AK 

UVPAR = ALUMUV*SINIJV*AK — 



C 


UVPUUB=UVBAP/UURAR 

H 820 




130 

c 


UV PUVR=UVPAP/( SORT (UUB A=> ) *SO FT (WRAP)) 

H 830 





c 


VVRUUR=SQRT CYVPAR) /SOPTCUU^A P) 

H 840 





c 


WRITE (6.40) HR.HI.UUBAR. WO OR ,U VP. AR , UVBUU8 , UV BUVR • VVRUUR 

H 85H 



- WRITE 16,40) HR.HI.UU3AP.WP AR . UVBAR , U BA PR , VBARR.OP 

IF CI.EQ.l) GO to 126 


135 



PU CI> = CUB APP-UL ) /(HT-VU 




— .... ... OUU ( I ) = C UUB AR-UUL ) / ( HI - V L ) 

0\MI) =CVPAPP-VL>/ (HT-YL) 

DV V ( I ) = ( V VP A R-VVL ) / C HI - Y L ) 

... DUVCI)=CUV«AR-UVL>/CHI-YL) 


14" 



PPPCI)=C OP-PPL) / CHI-YL ) 







l?g 

-VLSHI,, 




UUL=UUBAR 

VVlT.VV.PAR 

UVl^UVBAR 


_14f 



JJL =UR_APR „ 




VL=VRAPP 

nPL=OP 




1?9 

continue 

H 860 







.1 ?JL 

W D ITE (p>,60) UUOFTO 

H 9 1C 




15 f ___ . 



6K=12.0»A NH U 0*U 0 /^u 




- - - . .. .. . 

.HRrO,_Q-_P.HR . 

^0 200 1 = 1. IT 

hr=hr+ohr 

____ 

.15.5 



HI=HR»HU 










4 _o_)__hd , hi , duj < i > ,p v v m , puv (i ) , nil ci> .oven .pppcd 










Z£.!L 

CONTINUE 









return 

H 92P 
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END 
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